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for  Autonomous  Air  Vehicles 

Technical  Abstract 

Successful  operation  of  next-generation  unmanned  air  vehicles  will  demand  a  high  level 
of  autonomy.  Autonomous  low-level  operation  in  a  complex  environment  dictates  a  need  for  on¬ 
board,  robust,  reliable  and  efficient  trajectory  optimization.  In  this  report,  we  develop  and 
demonstrate  an  innovative  combination  of  traditional  analytical  and  numerical  solution 
procedures  to  produce  efficient,  robust  and  reliable  means  for  nonlinear  flight  path  optimization 
in  the  presence  of  time-varying  obstacles  and  threats.  The  trajectory  generation  problem  is  first 
formulated  as  an  optimization  problem  using  reduced-order  dynamics  that  result  from  the  natural 
time-scale  separation  that  exists  in  the  aircraft  dynamics.  Terrain  information  is  incorporated 
directly  into  the  formulation  of  the  reduced-order  dynamics,  which  significantly  reduces  the 
computational  load  and  leads  to  a  path  planning  solution  that  can  be  implemented  in  real-time. 
Various  cases  of  terrain,  pop-up  obstacles/threats,  and  targets  are  simulated.  A  representative 
optimal  trajectory  is  generated  with  in  a  high  fidelity  full-order  nonlinear  aircraft  dynamics  and 
compared  with  a  solution  obtained  from  a  reduced-order  optimization.  The  developed  algorithm 
is  flight  demonstrated  with  a  fixed-wing  unmanned  aircraft  test-bed  in  which  a  neural  network- 
based  adaptive  autopilot  is  integrated  with  the  on-line  trajectory  optimization  algorithm. 
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1.  Introduction 


1.1  Background 

High-flying  unmanned  reconnaissance  and  surveillance  systems  are  now  being 
used  extensively  in  the  United  States  military.  Current  development  programs  will  soon 
produce  demonstrations  of  next-generation  unmanned  flight  systems  that  are  designed  to 
perform  combat  missions.  In  practice,  these  vehicles  must  achieve  a  high  level  of 
autonomy  in  operations  to  be  successfully  deployed  in  large  numbers.  Their  use  in  first- 
strike  combat  operations  will  dictate  operations  in  densely  cluttered  environments  that 
include  unknown  obstacles  and  threats,  and  will  require  the  use  of  terrain  for  masking. 
The  demand  for  autonomy  of  operations  in  such  environments  dictates  the  need  for  an 
on-board  trajectory  optimization  capability. 

In  any  given  mission  scenario,  the  initial  conditions,  requirements,  objectives, 
vehicle  dynamics,  and  constraints  can  be  analyzed  at  varying  levels  of  detail  to  produce  a 
number  of  feasible  flight  plans.  However,  in  most  situations  conflicts  will  arise  between 
the  cited  factors  that  necessitate  trade-off  of  one  for  the  other.  Ad  hoc  methods  of 
solution  can  be  constructed  for  simple  missions  with  relatively  few  constraints,  but  such 
methods  are  quickly  overwhelmed  as  problem  complexity  grows.  In  such  cases  formal 
methods  for  trajectory  optimization  are  needed  to  identify  a  dynamically  feasible  solution 
that  is  “best”  in  some  sense. 
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Essentially  all  numerical  methods  for  solving  trajectory  optimization  problems 
incorporate  nonlinear  programming  (NLP)  methods.  Depending  on  how  NLPs  are 
employed,  methods  available  for  a  solution  of  trajectory  optimization  problems  generally 
fall  into  two  distinct  categories:  direct  and  indirect[l].  Indirect  methods  are  characterized 
by  explicitly  solving  the  optimality  conditions  stated  in  terms  of  the  adjoint  differential 
equations,  the  Pontryagin’s  maximum  principle,  and  associated  boundary  (transversality) 
conditions[2].  Using  the  calculus  of  variations,  the  necessary  conditions  are  derived  by 
setting  the  first  variation  of  the  Hamiltonian  function  to  zero.  The  indirect  approach 
usually  requires  the  solution  of  a  nonlinear  multi-point  boundary  value  problem.  An 
indirect  method  for  optimizing  a  function  of  n  variables  would  require  analytically 
computing  the  gradient  and  then  locating  a  set  of  variables  using  a  root-finding  algorithm 
such  that  the  gradient  is  zero.  In  contrast,  direct  methods  do  not  require  an  analytic 
expression  for  the  necessary  condition  and  typically  does  not  require  initial  guesses  for 
the  adjoint  variables.  Instead,  the  dynamic  state  and  control  variables  are  adjusted  to 
directly  optimize  the  objective  function[3].  All  direct  methods  introduce  some  parametric 
representation  for  the  control  variables  which  in  general  leads  to  larger  number  of 
variables  in  NLP. 

While  great  success  has  been  achieved  in  numerically  solving  complex  nonlinear 
trajectory  optimization  problems  using  direct  and  indirect  techniques,  many  short 
comings  still  exist.  Direct  optimization  schemes  tend  to  be  tolerant  of  a  poor  initial  guess 
but  are  in  general  slow  to  converge.  Indirect  methods,  in  contrast,  can  usually  achieve 
greater  accuracy  in  a  few  iterations  but  also  exhibit  much  greater  sensitivity  to  the  initial 
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guess.  Essentially  all  of  the  available  methods  in  both  categories  demand  computational 
resources  that  until  recently  were  simply  unavailable  in  a  flight  environment.  For 
complex  problem  fonnulations  with  full  vehicle  dynamics,  most  have  run  times  that 
preclude  real-time  mission  updates.  Almost  all  suffer  from  great  complexity  of 
implementation. 

For  these  reasons,  past  aircraft  trajectory  generation  studies  were  rooted  in 
classical  optimization  theory  and  the  calculus  of  variations.  A  great  variety  of  practical 
methods  for  optimizing  the  flight  of  aircraft  both  in  two  and  three  dimensions  were 
developed  in  the  1960s  and  70s.  Much  of  this  work  began  with  the  studies  performed  by 
Bryson  based  on  minimizing  time,  fuel  or  range  using  reduced  order  modeling 
methods[4,  5].  Subsequently,  the  use  of  singular  perturbation  methods,  and  multi-time 
scale  analysis  for  trajectory  optimization  were  suggested  by  Kelley [6,  7],  and  later 
extensively  developed  and  applied  by  Calise  [8-11],  Ardema  [12]  and  others.  In  the  work 
of  Calise,  the  emphasis  was  on  obtaining  analytic  or  near-analytic  feedback  solutions  that 
could  be  implemented  with  very  limited  computing  resources  in  real-time. 

In  this  research,  we  pay  attention  to  natural  time  scale  separation  that  exists  in 
aircraft  dynamics  and  attempt  to  exploit  two-time  scales  in  on-line  trajectory 
optimization.  The  overall  approach  follows  that  in  [13]  in  which  a  reduced-order  terrain¬ 
following  optimization  problem  is  formulated  for  the  nap-of-the-Earth  guidance  of 
helicopters  under  the  assumption  that  terrain-following  dynamics  are  effectively 
described  by  pseudo  3-D  dynamics.  Relying  upon  the  time-scale  separation  assumption, 
we  address  various  aspects  of  path-planning  such  as  moving  threats,  interior  point 


constraints,  and  pop-up  obstacles  besides  terrain-masking.  As  a  next  step,  we  implement 
the  resulting  optimal  trajectory  from  the  reduced-order  dynamics  in  a  high  fidelity  flight 
simulator  and  assess  the  resulting  trajectory  flown  in  the  high-fidelity  simulator.  To  be 
precise,  we  compare  the  trajectory  from  the  reduced-order  formulation  to  both  the 
trajectory  followed  in  the  high-fidelity  simulator  for  the  reduce-order  solution  and  the 
optimal  solution  for  full  order  6DOF  aircraft  dynamics.  In  that  procedure,  we  try  to 
reveal  various  aspects  that  should  be  taken  into  account  for  the  reduced-order  optimal 
trajectory  to  be  feasible  with  full-order  6  DOF  vehicle  dynamics.  The  report  includes 
flight-test  validation  to  demonstrate  on-line  re-planning  of  flight  trajectory  when  a 
simulated  pop-up  obstacle  appears. 

1.2  Contributions  of  the  PHASE  II  Efforts 

This  research  expands  on  the  work  done  by  Menon  and  Kim[13].  The  rationale  in 
this  research  is  to  resort  to  analytical  solutions  to  maximal  degree  possible  for  numerical 
efficiency.  With  this  philosophy,  it  was  decided  to  continue  investigating  using  a  reduced 
order  fonnulation  initiated  in  [13]  and  exploit  the  benefits  of  analytical  methods  of 
solving  the  optimal  path  planning  problem  before  it  became  necessary  to  use  any 
numerical  methods.  In  most  numerical  methods,  the  time  step  for  the  discretization  phase 
usually  must  be  very  small  to  avoid  loss  of  information.  This  leads  to  very  long  solving 
times.  Also,  the  graphical  methods  of  solving  this  type  of  problem  rarely  lead  to  an 
optimal  solution. 

An  underlying  theory  for  the  reduced-order  formulation  is  singular  perturbation. 
In  general,  singularly  perturbed  dynamic  systems  are  characterized  by  a  small  parameter 
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multiplying  the  derivatives  of  some  state  vector  components.  If  the  small  parameter  is  set 
to  zero,  the  order  of  the  dynamic  system  is  reduced.  The  solution  to  the  reduced  order 
problem  that  results  when  the  parasitic  small  parameter  is  set  to  zero  provides  the  leading 
term  in  a  perturbation  series.  Of  course,  the  reduced  order  solution  is  not  able  to  satisfy 
all  the  boundary  conditions  of  the  original  full  order  problem.  Depending  on  the  degree 
of  sub-optimality  of  the  reduced-order  solution  due  to  the  effect  of  neglected  dynamics, 
corrections  can  be  made  by  constructing  boundary  layer  transients  that  allow  rapid 
variation  of  the  fast  states  on  a  stretched  time  scale[14].  The  boundary  layer  solutions  are 
required  to  satisfy  the  boundary  conditions  violated  by  the  reduced  solution  and  to 
approach  the  reduced  solution  asymptotically.  Natural  time  scale  separation  that  exists  in 
aircraft  dynamics  allows  the  successful  application  of  singular  perturbation  methods  to 
problems  in  atmospheric  flight  path  optimization  [15].  The  reduced  and  boundary  layer 
problems  are  often  solvable  using  analytic  methods  alone.  This  has  made  a  unified 
analysis  of  3-D  high-performance  aircraft  path  optimization  possible  [16]. 

The  components  of  the  research  presented  here  include: 

•  Two  pseudo  3-D  formulations  that  include  wind  effects,  moving  targets, 
interior  point  constraints  in  the  form  of  waypoints,  and  moving  threats. 

•  Derivation  for  second  order  variation  conditions  for  each  formulation 
which  further  assists  finding  the  optimal  initial  conditions 

•  Expansion  of  pseudo  3-D  equations  of  motion  to  3-D  equations  of  motion 

•  Further  expansion  of  pseudo  3-D  equations  of  motion  to  3-D  ones  with 
addition  of  a  velocity  as  a  state.  This  constitutes  a  3-D  varying  velocity 
formulation. 
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•  Evaluations  for  the  ability  to  track  the  pseudo  3-D  optimal  trajectories  in 
a  high-fidelity  nonlinear  simulation  of  6  degree-of-freedom  (DOF)  high- 
performance  aircraft  operating  at  low  altitude  over  simulated  and  real 
terrains  with  both  stationary  and  moving  pop-up  threats 

•  Construction  for  the  6  DOF  full-order  optimal  trajectories  using  GESOP 
(Graphical  Environment  for  Simulation  and  Optimization),  the 
optimization  software  provided  by  University  of  Stuttgart. 

•  Comparison  and  analysis  of  the  full-order  and  pseudo  3-D  optimal 
trajectories 

•  Flight  demonstration  of  guided  flight  of  an  unmanned  aerial  vehicle 
(UAV)  while  employing  the  real-time  path  optimization  algorithm  by 
introducing  virtual  terrains  and  obstacles 


While  applying  corrections  for  the  optimal  trajectory  resulting  from  pseudo  3-D 
formulation  by  the  method  of  matched  asymptotic  expansion  as  in  [17]  were  proposed  in 
the  phase  II  proposal,  the  resulting  optimal  trajectories  from  pseudo  3-D  and  full  3-D 
formulations  did  not  exhibit  significant  differences,  both  of  which  resulted  in  trajectories 
that  exhibit  discrepancies  from  the  full-order  GESOP  optimal  solutions.  Various  issues 
for  comparing  pseudo  3-D  solutions  and  those  of  6  DOF  optimization  are  extensively 
addressed  in  [18]  the  work  of  which  was  supported  by  GST  for  this  Phase  II  effort  and 
performed  at  the  Georgia  Institute  of  Technology  under  the  direction  of  Prof.  Anthony  J. 
Calise.  The  analysis  using  the  singular  perturbation  method  for  full  6  DOF  dynamics, 
including  engine  modeling  and  aerodynamic  coefficients,  is  considerably  complex  and 
has  not  been  addressed  yet  in  the  literature  of  optimal  path  planning  for  aerial  vehicles 
that  generally  assumes  a  point  mass  model  for  optimization  problem.  This  task,  if 
undertaken,  would  require  extremely  complex  algebra  in  order  to  obtain  costate 
equations.  Therefore,  instead  of  attempting  to  carry  out  complex  singular  perturbation- 
based  6  DOF  optimization  analysis,  we  converted  our  attention  to  additional  aspects  of 
pseudo  3-D  formulation,  such  as  multi-vehicle  formulation  and  efficient  numerical 
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algorithms  for  the  case  of  increased  search  dimension  due  to  multi-vehicle  formulations. 
The  results  from  this  research  constitutes  the  thesis  of  Shannon  Twigg  [19].  So 
additional  contributions  from  the  Phase  II  efforts  are: 


•  Expansion  of  pseudo-3D  and  3D  equations  of  motion  to  handle 
cooperative  path  planning  for  multi-vehicles. 

•  Application  of  genetic  algorithms  (GA)  to  solve  multiple  initial 
conditions.  The  GA  employed  in  this  work  is  different  from  conventional  ones  as  in  [1] 
in  a  sense  that  GA  is  only  used  to  set  proper  initial  conditions.  With  a  reasonably  small 
number  of  variables  to  be  solved,  the  employed  algorithm  resulted  in  run-time  that  is 
comparable  to  the  variable  sweep  method  for  a  single  variable. 


12 


2.  Reduced-Order  Formulation 


In  this  section,  we  provide  two  pseudo  3-D  formulations  and  two  full  3-D 
formulations  extensively  studied  in  the  Phase  II  efforts.  For  ease  of  understanding,  we 
concentrate  our  presentation  on  those  of  pseudo  3-D  formulations.  Full  details  on 
expansion  of  the  methodology  to  full  3-D  formulations  and  that  of  multiple  vehicles  are 
referred  to  the  thesis[  1 9]  and  related  papers[20-24].  Figure  2.1  depicts  a  sample  terrain 
profile  with  the  X-Y-H  coordinate  system  and  a  local  xi-yi-zi  coordinate  system.  The 
moving  local  coordinate  system  has  its  origin  on  the  terrain  surface  at  a  current  x,  y 
position  with  the  xi-yi  plane  being  the  tangent  plane. 


x 


Figure  2.1:  Relationship  between  Inertial  Frame  and  Local  Tangent  Plane. 
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2.1  Local  Tangent  Plane  Equations  of  Motion 

The  local  tangent  plane  fonnulation  incorporates  the  constraint  that  the  vehicle 

flies  tangentially  to  the  local  terrain  directly  into  the  equations  of  motion  and  can  be 
written  as 


V cosy/  Kfxfy  sin^ 

A,  AxA2 


+  u{x,y) 


(2.1) 


y  = 


-  VAX  sin  y/ 

A 


+  v(x,y) 


(2.2) 


where  x  and  y  are  the  north  and  east  components,  respectively.  V  is  the  total  aircraft 
velocity  while  u  and  v  are  the  wind  velocities  in  the  x  and  y-directions,  respectively.  The 
heading  of  the  vehicle  is  represented  by  y/  —  the  heading  angle  measured  with  respect  to 
the  local  tangent  plane.  Also,/)  and  fy  are  the  partial  derivatives  of  the  terrain  profile.  A  / 
and  A  2  are  given  by 

A  =  t/1+/.2  (2-3) 

A  =  (l +  /.’+/, 2  (2.4) 


The  cost  function  for  this  problem  can  be  seen  in  the  following  equation. 

J  =  j()'  [0  -K)+  Kg(x,  y,  t)\lt 


(2.5) 
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In  this  equation,  the  combined  threat  and  terrain  function,  g(x,y,t),  is  given  as  a  function 
of  time  as  well  as  the  position  and  can  be  defined  as  follows. 

g(x,y,t)  =  f(x,y)  +  fr(x,y,t)  (2.6) 


Her Q,f(x,y)  is  the  function  for  the  terrain  profile  and  fr(x,y,  t)  is  the  function  denoting  the 
moving  threat.  The  weighting  parameter,  K,  can  vary  between  0  and  1  and  detennines 
the  relative  importance  of  time  and  terrain  masking/threat  avoidance  used  in  the 
optimization.  When  K  =  0,  the  equations  are  optimized  with  respect  to  time.  When  K  is 
set  to  1,  the  path  is  optimized  with  respect  to  the  threats  and  the  terrain.  The  Hamiltonian 
equation  can  then  be  given  as 


h  =  a4+ax 


V cos i//  VfJ  sin '// 

+ - - +  u 


A , 


AxA2 


+  A, 


-  VAX  sin  y/ 

aT~' 


+  v 


(2.7) 


In  this  expression,  Xx  and  Xy  are  the  costate  equations  and  A4  can  be  seen  in  the  following 
equation. 

AA=\-K  +  Kg{x,y,t)  (2.8) 


The  moving  threat  and  target  equations  of  motion  are,  respectively: 

xT  =  VT  cos  y/T 
yT  =  VT  sin  y/T 


(2.9) 


Xfg  V \g  cos  y/ jg 

vTg  =VTgsmy/Tg 


(2.10) 


In  each  expression,  it  is  assumed  that  the  respective  velocity  and  heading  angle  are 
known  at  all  times.  The  moving  target  then  results  in  a  new  boundary  condition. 
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¥(*/)  = 


(2.11) 


x(t)~xTg(t) 

y(t)-yTg(t) 


In  this  expression,  it  can  be  seen  that  *F(tf)  has  an  explicit  dependence  on  the  final  time  as 
a  consequence  of  the  fact  that  the  target  coordinates  are  assumed  to  satisfy  equation 
(2.10).  Therefore,  for  a  free  final  time,  the  Hamiltonian  equation  satisfies 

=  VTg  k  cos y/Tg  +  Av  sin  y/Tg  \=(  (2.12) 

t=tf 

Due  to  the  moving  threat,  the  Hamiltonian  equation,  (2.7),  is  explicitly  dependent  on 
time.  Given  this,  the  optimality  condition  for  a  solution  along  an  extremal  arc  shows  that 

H  =  —  =  Kgt  (2.13) 

ot 

where  gt  denotes  the  partial  derivative  of  the  penalty  function  with  respect  to  time. 
Assuming  that  the  threat  is  constant  when  expressed  in  a  coordinate  system  that  is 
attached  to  the  moving  threat,  then 

g(x,  y,  t )  =  g[x  -  xT  (t),  y-yT  (*)]  (2.14) 

with  the  threat  coordinates  satisfying  (2.9).  Thus 

H  =  -KVt (gx  cos if/T  +  gy  sin y/T)  (2.15) 

Because  the  final  time  is  free,  the  boundary  condition  for  this  expression  is  defined  in 

(2.12). 

The  optimality  condition  for  this  problem  is  defined  as 

Hv=  0  (2.16) 


H(tf)  =  -A1 


d y 

dt 
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Evaluating  this  expression  results  in  the  following  relationship 


^  =  K 


Vfjy 

A  A 


V  smy/ 

A 


A 

VAX  cos  y/ 


(2.17) 


Equation  (2.17)  can  then  be  substituted  into  the  Hamiltonian  equation,  (2.7),  to  determine 
equations  defining  the  two  costates,  A.x  and  Ay  as  follows. 


A  = 


-  ( A4  -  H)Ax  cos 


V 


Den 


A,  = 


(A4  -  H )A2  sin y/-(A4-H )ff  cosy/ 


Den 


where 


Den  =  VAX  +  Afu  cos  y/  +  fxfyv  cos  y/  -  A2v  sin  y/ 


(2.18) 


(2.19) 


(2.20) 


These  new  expressions  for  the  costates  can  then  be  inserted  into  (2.12)  to  result  in  a  new 
boundary  condition  for  the  Hamiltonian  at  the  final  time. 


VTg  A4  (ax  cos  y/  cos  y/Tg  +  fxfy  cos  y/  sin  y/Tg  -  A2  sin  y/  sin  y/Tg ) 
VT  [ax  cos  yr  cos  y/Tg  +  fxf  cos  y/  sin  y/Tg  -  A2  sin  y/  sin  y/Tg )  -  Den 


(2.21) 


Differential  equations  for  the  costates  can  be  found  using 

K  =  ~H* 

K  =  ~Hy 


(2.22) 


This  yields 


A  =  ~Kgx-A 


D2  cos y/  +  D2  sin  y/  +  Dxux 

D 


-A 


D4siny/  +  Dxvx 

A 


(2.23) 
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(2.24) 


Av=-Kgv-Ax 


D5  cos  y  +  A  sin  y/  +  D{u^ 

A 


-A, 


D-j  sin  y/  +  D,v 


1  vy 


D, 


where 


Dl  =  A?  4 

(2.25) 

A  =  -vfxf„A 

(2.26) 

VA\A\B2  -  VA;4fyfxx  -  VA?BJxfy 

(2.27) 

D4=VA1%-VAAfxf« 

(2.28) 

A  =  -VfxfxyA\ 

(2.29) 

VAAB3 -V44fJ\v-lA'BJ\fv 

(2.30) 

A  =  VA,4B4  -  VAf  A2fxfxy 

(2.31) 

4  =  fj „  +  fyfv 

(2.32) 

A  =  fjv  +  fyfxx 

(2.33) 

A  =  f  x  f  yy  +  f  y  f  xy 

(2.34) 

A  =  fxfxy  +  fyfyy 

(2.35) 

Next,  the  time  derivative  of  either  equation  (2.18)  or  (2.19)  is  taken  and  set  equal  to  its 
counterpart  in  equation  (2.23)  or  (2.24).  This  expression  can  then  be  solved  for  the 
derivative  of  the  heading  angle  such  that 
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where 


¥  = 


Tx  +  T2u  +  T3v  +  T4ux  +  T5u  +  T6vx  +  T7v 


7  y 


(2.36) 


7)  =  -KVSX  +  V(A4  -  H)S2  (2.37) 

T2  =-KS3  +(A4  -H)s4  (2.38) 

T3=KS5+(A4~h)S6  (2.39) 

T4=(A4-H)S7  (2.40) 

Ts  =(A4  ~H)Ss  (2.41) 

T6  =  (A4  -  H )S9  (2.42) 

T7={A4-H)Sxo  (2.43) 

T%  =  (A4  ~  H )a]  A\  (2.44) 

Sx  =  A?A2[A2gx  sin y/  +  (gy  +  fx2gy  - gJJy )cosi//]  (2.45) 

S2  =  (fjyfxx  -  A  |2  f y  f xy  )sin  ¥  +  Afyftx  COS  <//  (2.46) 

53  =  43^2  sin  ^  +  (gy  +  fx2gy  -  gJJy  )cos ^]cos ^  (2.47) 

S4=AlA2(A2fxfxy  -///,/„  +/,/„)■  cos>  + 

/  ~  o  o  \  (2.48) 

+  A  [A  fxf„  -  A  f y  f xy  +  2  fxfy  /„  )SI11  ¥  COS  I// 

=  A 4  U2 g  v  -  2SjJy  )sin  igcosig+  ^  ^ 

+  AAlgx  sin2  ¥  +  AAi{gxfxfy  ~  42/v^v/v)cos2 1// 
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(2.50) 


$6  =  AA2(fyfxy  -flfyfv  +  47,/Jcos2  V  + 

+  A  (2  fjyfxy  +  •  .4  )sin  V  COS  I// 

4  =  7 4  (4  sin  y/  -  fxfy  cos  yr )cos  yr  (2.51) 

S8  =  A\A2  cos2  yr  (2.52) 

4  =  AAl(2fJy  sin^cos y/ -  A2)+  AxA2{a22  -/t2/;)cos2  y/  (2.53) 

sw  =  AlAi(fJy  cos y/-  A2  sin^)cos^  (2.54) 


This  solution  consists  of  four  differential  equations  —  x,  y,  H  and  y/  —  and  requires  two 
initial  conditions  to  be  found  —  H  and  y/.  The  final  value  of  the  Hamiltonian  is  known, 
via  equation  (2.21).  The  solution  is  reached  when  the  final  values  of  the  Hamiltonian  and 
position  are  met  and  the  cost  is  minimized.  When  there  are  no  moving  threats,  the 
Hamiltonian  is  constant  in  value  -  so  there  are  only  three  differential  equations  -  and  the 
final  value  is  still  known.  When  there  is  no  moving  target,  the  final  value  of  the 
Hamiltonian  is  zero. 


2.1.1  Legendre-Clebsch  Necessary  Condition 

The  Hamiltonian  equation  for  the  local  tangent  plane  equations  of  motion  is 


H  =  A,+K 


V  cos  y/  +  VfJySmy/ 


A 


A  A 


+  A. 


-  FA,  sin  y/ 

Z 


(2.55) 


and  the  algebraic  equations  for  the  costates  are 
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(2.56) 


4 


~(A4  -  H)Ax  cos  y/ 

V 

(A4  -  H)(a2  siny/  -  fjy  cos  y/ ) 


VA\ 


Since  there  is  only  one  control,  the  second  partial  derivative  of  the  Hamiltonian  with 
respect  to  the  heading  angle  is  a  scalar  value  and  is  represented  by 


-  V  cos  y/ 

VfJ\,smyy~ 

+  4 

VAl  sin  y/ 

A 

A\A2 

A 

Huu  =  4 

Substituting  in  the  equations  for  the  optimal  costates  will  result  in 

huu={a4-h)>  0 


(2.57) 


(2.58) 


This  condition  must  always  be  satisfied. 


2.1.2  Weierstrass  Test 

The  variational  Hamiltonian  can  be  found  by  substituting  the  costate  equations 
from  (2.56)  for  the  optimal  path  into  the  Hamiltonian  equation  from  (2.55)  evaluated  for 
any  path.  This  yields 


h(w)  =  a4  + 


~{A4  -H)Ax  cos y/° 

V cosy/  |  Vfxfy  sin  y/ 

V 

A\  Ax  A2 

(A4  -  H )(a2  sin  y/  -  fjy  cos y/° ) 

-  VAX  sin  y/ 

VA 

A 

(2.59) 


This  can  be  simplified  to 

H(V/)  =  (A4  -7/)[l-cos(^°  -y)]>0  (2.60) 

which  will  always  be  satisfied  if  equation  (2.58)  is  satisfied. 
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2.2  Simplified  Equations  of  Motion 

The  equations  of  motion  used  in  the  simplified  formulation  are: 


x  -  V cosy /  +  u(x,y ) 

(2.61) 

y  =  V  sin^  +  v(x,y) 

(2.62) 

These  equations  are  written  in  the  local  level  plane  and  neglect  the  effects  of  the  terrain 
slope.  The  cost  equation  for  this  case  is  the  same  as  earlier  and  can  be  found  in  equation 
(2.5).  The  corresponding  Hamiltonian  equation  is  therefore 


H  =  A4  +  Xx  [V  cos  y/  +  u]  +  Xy  [V  sin  y/  +  v] 


(2.63) 


The  equations  governing  the  moving  target  and  moving  threat  can  be  seen  above  in 
equations  (2.9)  and  (2.10).  Evaluating  the  optimality  condition  stated  in  equation  (2.16) 
for  this  fonnulation  results  in  the  expression 


K  =  K 


sin^ 

cos^ 


(2.64) 


Substituting  this  into  the  Hamiltonian  equation  results  in  the  following  costate  equations 


-  (A4  -ff)cosy/ 

V  +  u  cos y/  +  vsiny/ 


(2.65) 


X„ 


~(Aa  -H)siny/ 

V  +  a  cos y/  +  vsin  y/ 


(2.66) 


Therefore,  the  Hamiltonian  evaluated  at  the  final  time  will  be 


VTgA4  cos(^-^rg) 

VTg  cos(^  -  y/Tg )  -  (V  +  u  cos  y/  +  v  sin  y/j 


(2.67) 
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The  costate  differential  equations  can  then  be  found  to  be 


K  =  ~H ,  =  -Kg, 
K=~Hy=-Kgy 


AxUx~AyVx 

(2.68) 

-^Vy 

(2.69) 

As  before,  the  time  derivative  of  (2.65)  or  (2.66)  is  found  and  equated  to  either  (2.68)  or 
(2.69).  This  expression  can  then  be  rearranged  to  result  in  the  following  heading 
differential  equation. 


¥  = 


R\  +  +  RiV  +  R4  (»  v  -  Vy)+R5Uy  +  Vx 

Ri 


(2.70) 


with 


=  KV(gv  cos  y  -  gx  sin  y ) 

(2.71) 

^2  =  K  (g ,,  cos  y  -  gx  sin  y  )cos  y 

(2.72) 

R3  =  K[gy  cos  y  -  gx  sin  y  )sin  y 

(2.73) 

R4  =  (A4  -H)siny  cos  y 

(2.74) 

R5  =  ~(A4  -  H) cos2  y 

(2.75) 

R(>  =  (A4  -H) sin2  y 

(2.76) 

I 

II 

oT 

(2.77) 

Again,  the  inclusion  of  a  moving  target  and  moving  threat  results  in  a  system  of  four 
differential  equations  with  two  initial  parameters  to  be  found. 
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2.2.1  Legendre-Clebsch  Necessary  Condition 

The  Hamiltonian  equation  for  the  simplified  equations  of  motion  is 

H  =  A4  +  Ax  \V  cos  y/  ]  +  Xy  [V  sin  y/  ] 


and  the  costate  equations  are 


K 

K 


-{A4  -H) cost// 

V 

-  (A4  -H)sini// 

V 


(2.78) 


(2.79) 


The  partial  derivative  of  the  Hamiltonian  equation  with  respect  to  the  heading  angle  is 


h,.=(a,-h)>0 


(2.80) 


whish  must  always  be  satisfied. 

2.2.2  Weierstrass  Test 

Using  the  Hamiltonian  equation  in  (2.78)  and  the  costate  equations  in  (2.79),  the 
variational  Hamiltonian  can  be  written  as 


H  =  A4  + 


~(A4  -H) cost//' 
V 


[i V  cos^]- 


~(Aa  -H)sini//C 
V 


[i V  sin^] 


(2.81) 


This  equation  can  be  reduced  to 

H{y/)  =  (A4  -H)[\-  cos(y/°  -^)]>0  (2.82) 

which  will  always  be  satisfied  if  equation  (2.80)  is  satisfied. 


2.3  Interior  point  constraints  (way  points) 
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One  useful  problem  variation  that  can  also  be  investigated  is  the  implementation 
of  interior  point  constraints  during  the  flight.  For  this  problem,  the  constraints  will  be 
implemented  in  the  form  of  waypoints  where  a  specific  position  is  required  in  the  middle 
of  the  flight.  There  can  be  n-number  of  waypoints  during  this  flight,  such  that  each 
waypoint  -  with  a  given  x  and  y  position  -  is  reached  at  an  unspecified  time,  th  in  a 
specified  order  before  ending  at  the  specified  final  position. 

In  this  type  of  problem,  there  are  certain  constraints  on  the  costates  and 
Hamiltonian  that  must  be  fulfilled  at  the  interior  points.  They  include 

4(h+)  =  )  +  uu 

Ay(t;)  =  Ay(t;)  +  v2i  (2.83) 

H{tt)  =  H{t:)  =  0 


This  implies  that  the  value  of  each  of  the  costates  will  jump  at  each  waypoint  while  the 
Hamiltonian  will  remain  constant.  Because  of  that,  the  heading  angle  will  also  jump  at 
each  waypoint.  This  will  result  in  a  trajectory  such  as  that  seen  in  Figure  2.2.  In  this 


1.3 

1.2 

1.1 

1 

0.9 


0.7  ■ 

0.6  - 
0.5  ■ 

0.4  : _ 

0.3  I - 1 - 1 - 1 - 1 - 1 - 

0  2  4  6  8  10  12 

t  (sec) 


Figure  2.2:  Solution  with  one  waypoint. 
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example,  flight  over  a  flat  plane  with  one  waypoint  is  considered.  Because  all  the  terrain 
partial  derivatives  are  zero,  the  local  tangent  plane  equations  of  motion  will  reduce  to  the 
simplified  equations  of  motion.  In  addition,  both  the  costates  will  be  at  a  constant  value 
at  all  times  with  a  jump  at  the  time  of  reaching  the  waypoint.  This  results  in  the  heading 
angle  also  being  a  constant  value  with  a  jump  at  the  waypoint.  Having  a  jump  in  the 
heading  angle  will  create  an  optimal  path  that  is  not  flyable.  Therefore,  the  equations  of 
motion  for  this  section  will  be  modified  to  ensure  a  smooth  trajectory. 

The  equations  of  motion  for  this  section  will  include  the  equations  used  earlier 
with  the  addition  of  i//as  an  additional  state.  This  results  in  equations  of  motion  of 


.  V cosy/  Vfxf  sin y/ 

X  = - —  H - 1 - 

A\  AxA2 
-  VA,  sin  i// 

y  =  — 1 - 

a2 

y/  =  u 


(2.84) 


for  the  local  tangent  plane  or 

x  =  V  cosy/ 

y  =  Vsmy/  (2.85) 

y/  =  u 


for  the  simplified  equations  of  motion.  In  these  equations,  u  designates  the  control 
variable  for  the  system.  The  new  cost  equation  is 

J  =  \‘f  {(l - K  +  Kf)+Wu2}dt  =  \'f  {a4  +  Wu2)dt  (2.86) 


which  is  the  same  as  above  with  the  inclusion  of  the  control  in  the  cost  multiplied  by  a 
weighing  factor. 
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The  Hamiltonian  for  this  problem  is  now 


H  —  +  JVu  +  Axx  +  Avy  +  Xjf u 

For  both  fonnulations,  evaluating  the  optimality  equation  results  in 

Hu  =  0  =  2Wu  + 

This  yields  the  following  equation  for  the  control 

• 

u  =  w  = - 

2  W 


(2.87) 


(2.88) 


(2.89) 


At  each  interior  point,  the  following  conditions  on  the  costates  and  the  Hamiltonian  must 
be  met. 


(K  )  -  )  +  yi/ 

K  (*! )  =  ky  (/, ) + o2i 

H(t;)  =  H{t:)  =  0 


(2.90) 


Using  these  conditions  from  (2.90)  as  well  as  equations  (2.87)  and  (2.86),  an  independent 
equation  for  Xy  is  found  as  follows. 


A 


y 


Wu 1  -  A4  -  Axx 

y 


(2.91) 


Next  the  differential  equations  for  the  other  two  costates  can  be  determined  using 


K  =  ~Hx 

K  =  ~H r 


(2.92) 


This  yields 
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(2.93) 


Kfx-Ax 

D2  cosy/  +  D3  sin^ 

-A 

D4  sin  y/ 

A 

j 

D 

1 

j 

Vsiny/  Vfxfy  cosy/ 

+  Ay 

VAl  cosy/ 

lx 

A 

A{A2 

a2 

for  the  local  tangent  plane  equations  of  motion  with  Dj  through  D7  defined  in  equations 
(2.25)  -  (2.31).  For  the  simplified  equations  of  motion,  the  costate  differential  equations 
are 


K  =  -Kf* 

X  =  AXV  sin  y/  -  AyV  cos  y/ 


(2.94) 


This  results  in  a  system  of  six  differential  equations.  The  initial  conditions  for  the  two 
costates  -  Ax  and  Av-  must  be  found  as  well  as  Uj  for  each  waypoint. 


2.3.1  Legendre-Clebsch  Necessary  Condition 

The  Hamiltonian  equation  is  stated  in  equation  (2.87).  The  second  partial 

derivative  of  it  then 


Huu  =2W>0 

which  means  that  this  condition  is  always  satisfied. 


(2.95) 


2.3.2  Weierstrass  Test 

Using  the  Hamiltonian  equations  in  (2.87)  and  the  algebraic  equations  for  Ay  and 
Av  found  in  (2.91)  and  (2.89),  the  variational  Hamiltonian  can  be  found  to  be 

H(u)  =  w(u-u°)2  >  0  (2.96) 
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This  shows  that  the  Weierstrass  test  is  always  satisfied. 


2.4  Expansion  to  3-D  formulations 


The  3D  equations  of  motion  used  are 

x  =  V  cosy  cos  y/ 
y  =  V  cos  y  sin  y/ 


z  =  V  sin  y 


(4.1) 

(4.2) 

(4.3) 


Here,  V  is  the  vehicle  velocity,  y  represents  the  flight  path  angle  and  y/  is  the  heading 
angle.  The  cost  equation  used  for  this  problem  is 

J  =  J',{C,+C2}rf/  (4.4) 

Ct=\-K  +  Kf(x,y)  (4.5) 

C2=W[z-(f(x,y)  +  hJ  (4.6) 

This  cost  equation  has  two  distinct  parts.  The  first,  and  dominant  part,  is  Ci  shown  in 
(4.5).  This  part  controls  the  importance  of  minimizing  terrain  masking  versus 
minimizing  flight  time.  The  second  part  is  C2  as  seen  in  (4.6).  Here,  hc,  the  ground 
clearance,  is  a  constant  provided  by  the  operator  and  represents  the  desired  flight  height 
above  the  terrain.  This  part  is  used  to  keep  the  flight  path  near  the  desired  ground 
clearance  throughout  the  flight.  W  is  a  weighing  parameter  supplied  by  the  user. 
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In  case  the  velocity  of  the  vehicle  is  included  as  an  additional  state,  the  new 
equations  of  motion  will  now  consist  of  the  equations  in  (4. 1-3)  and  also 

•  T-D 

V  = - gsiny  (4.27) 

m 

In  this  equation,  m  is  the  mass  of  the  vehicle,  g  is  gravity,  and  T  is  the  thrust  of  the 
vehicle,  which,  for  optimal  results  should  be  held  constant  at  its  maximum  value.  The 
optimization  procedure  for  this  case  is  detailed  in  [19]  and  thus  omitted  here. 
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3.  Six  DOF  Optimization  by  GESOP 


3.1  Modeling 

The  6  DOF  full-order  dynamics  for  GESOP  optimization  are  constructed  by 
incorporating  engine  dynamics  and  propeller  rotation  that  are  built  into  the  f-wing- 
simulator  (Simulator  programmed  and  used  by  the  Flight  Mechanics  and  Controls  Group 
at  Georgia  Tech.).  This  f-wing  simulator  is  utilized  to  evaluate  the  feasibility  of  the 
pseudo  3-D  optimal  trajectory  as  well.  The  6  DOF  simulation  dynamics  comprise  14 
states  and  described  by 


r=TTBV  (1.7) 

mV  =  -ma>xV  +  F  (1.8) 

a)  =  -Il(a>xIa))  +  IlM  (1.9) 

«  =  (i.io) 


where  r  eR3  is  the  position  vector,  V  e  R3  is  the  velocity  vector,  cb  e  R}  is  the  turn 
rate,  and  q  e  R4  is  the  quaternion  vector,  TB  e  R  ''  and  Clq  e  7? 4x4  are  the  directional 
cosine  matrix  and  the  quaternion  matrix,  m  is  the  mass  of  the  vehicle,  and  /  is  the 
moment  of  inertia.  The  force  term  F  and  the  moment  term  M  are  detennined  by 
aerodynamic  coefficients  and  engine  characteristics.  They  are  based  on  the  Pioneer  UAV 


31 


with  some  changes  in  the  physical  dimensions.  The  calculation  of  the  aerodynamic 


coefficients  for  f-wing  and  the  6DoF-model  is  based  on  the  coefficient  tables  in  [25],  For 
the  GESOP-6DoF-Model  these  coefficient  tables  needed  to  be  curve  fitted  using  splines 
to  provide  smooth  functions  to  the  optimizer  without  a  significant  loss  of  performance. 
Table  3.1  gives  an  overview  over  the  basic  model  characteristics  used  for  the  6  DoF- 
Optimizations. 


Si-system 

Non-SI-system 

Wing  area 

1.8910  m2 

20.3546  ft2 

Wing  width 

4.3160  m 

14.1600  ft 

Chord  width 

0.4382  m 

1.4375  ft 

Air  density 

1.2250  kg/m3 

2.3770  slugs/ft3 

Mass 

65.4  kg 

4.4813  slugs 

Earth  acceleration  g 

9.8067  m/s2 

32.174  ft/s2 

Moments  of  inertia 

I XX 

7.0838  kg  m2 

5.2248  slugs  ft2 

Iyy 

13.6422  kg  m2 

10.062  slugs  ft2 

Izz 

16.7213  kg  m2 

12.333  slugs  ft2 

Ixy 

0.0  kg  m2 

0.0  slugs  ft2 

Ixz 

-0.9965  kg  m2 

-0.735  slugs  ft2 

lyz 

0.0  kg  m2 

0.0  slugs  ft2 

32 


Table  3.1:  Model  characteristics 


3.2  Issues  for  comparison  with  3-D  results 

Recall  that  the  pseudo  3-D  formulations  in  Sections  2.1  and  2.2  incorporate  the 
vertical  motion  into  the  equations  of  motion  so  that  the  resulting  trajectory  maintains  the 
same  altitude  over  the  terrain.  The  3-D  formulation  in  Section  2.4  or  6  DOF  formulations 
in  Section  3.1,  raise  the  question  of  how  to  distinguish  hills  to  prevent  the  optimizer  from 
leading  to  a  solution  that  flies  through  hills.  One  possibility  is  put  a  high  penalty  on 
changing  the  altitude  over  the  ground  so  that  the  optimal  path  lays  somewhere  close  to 
the  constant  altitude  over  the  ground  layer,  which  is  precisely  what  has  been  done,  as  in 
the  term  C4  in  (4.6),  in  the  3-D  formulation  in  Section  2.4.  In  case  of  the  6  DOF  full 
aircraft  dynamics,  since  the  aircraft  dynamics  involve  more  states  and  control  variables 
neglected  in  the  reduced-order  formulations,  ensuring  non-contact  between  the  airplane 
and  the  terrain  should  be  accomplished  by  properly  maneuvering  the  airplane,  which  in 
general  requires  controlling  the  aircraft  by  control  effectors  and  throttle.  This  point 
indicates  that  direct  comparison  between  optimal  trajectories  from  the  reduced-order  and 
the  full  6  DOF  aircraft  dynamics  is  not  straightforward.  For  example,  even  if  the  same 
cost  functional  is  used,  whereas  the  pseudo  3-D  formulation  leads  to  a  trajectory  of 
constant  velocity,  the  6  DOF  optimizer  cannot  be  forced  to  fly  at  a  constant  velocity  and 
the  cost  integral  over  time  will  have  little  meaning  for  direct  comparison. 

This  indicates  that  for  a  full  6  DOF  optimal  trajectory  to  be  comparable  to  the 
reduced-order  solution,  we  need  to  address  the  following  issues: 

•  How  to  ensure  for  the  flight  path  in  the  6  DOF  optimizer  to  be  close  to  the 
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optimal  trajectory  of  the  reduced-order  formulation  so  that  the  cost 
functional  in  (2.5)  makes  sense  in  both  cases 

•  How  to  set  boundary  conditions  for  extra  states  in  addition  to  those  states 
resulting  from  the  reduced-order  dynamics. 

•  How  to  transform  constraints  on  the  control  effectors  and  the  throttle  into 
equivalent  constraints  on  those  states  that  are  treated  as  controls  in  the 
reduce-order  fonnulation. 

Case  studies  for  these  complexities  and  some  related  simulations  are  presented  in  [18]. 
For  example,  an  experiment  with  the  terrain  altitude  as  a  path  constraint  in  GESOP  full 
order  optimization  required  a  large  number  of  grid  points  for  real  terrain  to  cover  all 
terrain  peaks.  Also,  direct  optimization  of  6  DOF  model  produced  heavy  peaks  and  jumps 
in  control  signals,  which  are  completely  neglected  in  the  reduced-order  formulation. 
Therefore,  to  circumvent  those  complexities  arising  due  to  aircraft  dynamics  such  as 
control  peaks  and  engine  characteristics,  additional  cost  terms  are  included  in  the  6  DOF 
optimization  procedure[18].  The  goal  for  these  additional  costs  is  to  realize  the  constant 
altitude  flight  path  in  the  full  aircraft  dynamics  as  close  as  possible  to  those  in  the 
reduced-order  formulation.  The  additional  terms  in  GESOP  optimization  includes 
avoiding  peaks  in  controls,  altitude  control  by  terrain-following,  and  inhibiting 
trajectory-terrain  intersection.  The  additional  cost  induced  by  these  terms  for  terrain- 
masking  is  calculated  and  compared  to  the  original  cost  in  a  terrain-following  scenario  as 
well  in  [18]. 
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4.  Numerical  Procedure  for  Optimal 

Trajectories 


4.1  Optimization  Procedure  for  the  reduced-order  formulation 

In  case  of  reduced-order  formulation,  a  solution  that  satisfies  a  set  of  necessary 
conditions  for  optimality  is  found  by  seeking  proper  initial  conditions.  Two  different 
methods  were  utilized  to  find  these  initial  conditions,  depending  on  the  number 
necessary.  When  only  one  value  was  needed,  a  variable  step  sweep  was  employed  to  find 
it.  Otherwise  the  GA  was  used[26-28]. 

To  begin  the  genetic  algorithm,  a  set  of  48  chromosomes  was  initialized 
representing  different  sets  of  initial  conditions  to  test.  Each  initial  value  in  the 
chromosomes  was  represented  by  digits  with  five  decimal  places  included.  In  addition 
the  costates,  flight  path  angles  and  Hamiltonian  values  included  an  extra  digit  to  indicate 
a  positive  or  negative  value.  After  the  chromosomes  were  initialized,  they  were  each 
tested  to  detennine  their  relative  costs.  To  accomplish  this,  the  current  chromosome 
being  tested  was  broken  into  its  respective  initial  conditions,  which  were  then  used  in  the 
differential  equations.  The  cost,  J,  was  found  for  the  run  as  well  as  the  distance  from  the 
final  position  of  the  run  to  the  final  target  position.  The  sum  of  these  two  values  was 
used  as  the  total  cost  for  the  chromosome. 

After  each  chromosome  was  tested  and  a  total  cost  assigned,  the  chromosomes 
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were  ranked  from  last  to  first  based  on  a  tournament  procedure.  Two  chromosomes 
would  be  randomly  chosen  to  compete  and  the  one  with  the  higher  cost  was  placed  in  the 
next  position  on  the  list  while  the  one  with  the  lower  cost  was  returned  to  the  available 
set  of  chromosomes  to  be  tested.  This  process  was  continued  until  all  the  chromosomes 
were  ranked.  The  top  24  chromosomes  were  then  kept  to  begin  the  next  generation. 

In  all  subsequent  generations,  the  24  available  chromosomes  were  combined  to 
create  24  new  chromosomes  to  complete  the  population  of  48.  Here,  chromosomes  1  and 
2  would  be  combined  to  create  two  new  chromosomes,  and  then  chromosomes  3  and  4 
would  be  combined  to  create  two  new  chromosomes  and  so  on  until  all  the  chromosomes 
were  mixed.  This  was  accomplished  by  first  mixing  the  individual  segments  of  the 
chromosomes  so  that  each  of  the  new  chromosomes  had  some  segments  from  each 
parent,  where  a  segment  consisted  of  the  digits  for  each  initial  condition  needed.  Next  a 
mutation  was  introduced  into  the  new  chromosomes  such  that  up  to  about  a  third  of  the 
digits  could  be  changed.  The  number  of  digits  changed,  which  digits  were  changed,  and 
their  new  values  were  all  determined  randomly.  After  all  the  new  chromosomes  were 
created,  the  cost  assignment  and  tournament  were  repeated  as  before.  This  process  was 
repeated  until  it  converged  on  a  solution. 

The  differential  equations  were  solved  using  a  standard  fourth  order  Runge-Kutta 
method.  In  addition,  a  variable  time  step  was  implemented  to  decrease  the  time  needed 
to  numerically  solve  the  set  of  differential  equations.  For  most  of  the  flight,  the  time  step 
was  0.1  seconds;  however,  when  the  distance  to  the  target  final  position  was  close 
enough,  the  time  step  was  decreased  to  0.01  seconds. 

Another  condition  was  added  to  the  differential  equation  solver  to  help  decrease 
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the  solving  time  of  the  genetic  algorithms.  For  each  formulation,  inequalities  were 


derived  that  had  to  be  satisfied  at  all  times  in  order  for  the  Legendre-Clebsch  necessary 
condition  and  the  Weierstrass  test  to  be  satisfied.  These  inequalities  were  then  tested  at 
each  time  step.  If  either  was  violated,  then  the  current  run  was  ended  at  that  point.  This 
decreased  the  solving  time  significantly,  but  was  even  more  useful  in  ensuring  the 
convergence  to  a  strong  local  minimum. 

Table  4.1  contains  some  average  run  times  for  the  different  formulations  using  1.8 
GHz  Pentium  IV  with  1  GB  RAM.  This  table  contains  the  number  of  initial  conditions  to 
be  solved  for,  the  run  time  length  for  each  problem,  then  the  average  time  it  took  to  solve 
the  problem.  It  can  be  seen  that  the  single  vehicle  formulations  were  all  generally  solved 
in  less  than  a  minute  for  a  case  with  a  run  time  of  10  seconds.  The  multiple  vehicle 
formulations  took  longer  to  solve  because  those  cases  tended  to  have  a  large  number  of 
local  minima. 


Fonnulation 

Number  of  Variables 

Run  time 

Time  to  Solve 

Pseudo  3D 

1 

10  sec 

25  sec 

3D 

2 

10  sec 

40  sec 

Varying  Velocity 

3 

10  sec 

65  sec 

2-vehicle 

3 

25  sec 

6  min 

3 -vehicle 

5 

25  sec 

15  min 
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4-vehicle 

7 

25  sec 

60  min 

2-vehicle  3D 

5 

10  sec 

10  min 

Table  4.1:  Time  to  Solve  for  Reduced-order  Formulations 


4.2  Optimization  Procedure  for  GESOP 

GESOP  provides  several  different  optimizers.  The  two  main  types  are  methods 
either  using  a  Direct  Collocation  Method  or  a  Direct  Multiple  Shooting  Method.  All 
optimizations  in  this  work  were  done  with  Multiple  Shooting  Method  called  SNOPT 
(Sparse  Nonlinear  OPTimizer)[29]  . 

The  maximum  number  of  iterations  was  changed  from  the  default  80  to  much 
higher  values  for  the  more  complex  iterations.  Values  around  1000  were  used  to  make 
sure  that  the  optimizer  was  able  to  find  the  optimal  solution  without  any  restarts.  Initially 
restarts  were  used  to  obtain  a  stable  solution,  but  later  it  was  found  that  the  majority  of 
cases  did  not  require  a  restart. 

The  values  for  “Real  Workspace  Size”  and  “Integer  Workspace  Size”  had  to  be 
increased  up  to  1,000,000  and  300,000  for  the  most  complex  optimizations,  to  inhibit  the 
raise  of  SNOPT  error  code  20  which  indicates  a  too  small  workspace.  The  variables 
“Opt.  Tolerance”  and  “Const.  Tolerance”  were  left  atl.OxlO-6.  This  causes  relatively 
high  iteration  times  but  on  the  other  side  a  very  accurate  solution.  The  change  of  the 
solution  on  the  last  30%-50%  of  all  iteration  steps  often  consists  of  a  very  low  and  almost 
negligible  change  over  a  good  part  of  the  last  iteration  steps.  For  the  Major  Grid  the 
number  of  points  was  increased  to  50  for  most  cases.  Especially  the  optimizations  for 
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flights  over  “Real  Terrain”  needed  such  high  values  to  converge.  The  Constraints  Grid 
was  set  to  a  single  point  since  for  the  final  cases  no  constraints  were  applied.  The 
Controls  Reference  Grid  was  set  to  a  relatively  low  value  of  9  points  in  most  cases.  For 
the  simulation  of  the  complex  optimizations  the  “Number  of  Points”  for  the  output 
spacing  had  to  be  increased  to  201. 

The  average  run  time  needed  to  converge  is  mostly  dependent  on  the  accuracy 
required.  With  a  very  high  accuracy  setting  the  optimal  solution  is  found  in  most  cases 
after  about  45-90  minutes  CPU-time  (1.8Ghz  Pentium  IV).  With  a  lower  accuracy  setting 
the  optimal  solution  can  be  found  in  about  30-45  minutes  with  just  very  little  differences 
compared  to  the  high  accuracy  setting.  A  good  result  which  approximately  follows  the 
optimal  solution  is  available  already  after  a  couple  of  minutes.  The  results  of  each  step 
can  be  observed  already  during  the  optimization  process  within  GESOP.  In  general 
optimizations  going  over  real  terrain  need  some  more  time  compared  to  mathematical 
terrains  due  to  the  interpolation  routine. 
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5.  Numerical  Results 


5.1  Terrain  Data 

Two  general  types  of  terrain  models  are  used  for  the  results  throughout  this 
report.  The  first  is  a  generic  terrain  model  used  for  the  initial  testing  of  the  equations. 
This  consists  of  variations  of  a  flat  plane  with  one  or  more  constructed  hills.  The  second 
consists  of  actual  terrain  data  for  a  larger  area.  This  allows  the  opportunity  for  the 
various  equations  to  be  tested  in  a  more  realistic  manner. 

A  sample  terrain  of  the  generic  model  is  shown  in  Figure  5.1.  In  this  case  a 
mostly  flat  plane  with  a  single  hill  is  used.  This  hill  in  this  terrain  is  formulated  using  the 
exponential  function 

f  =  Aerl/b  (5.1) 
where  A  is  the  amplitude,  b  is  a  scaling  factor  to  adjust  the  width  and  r  is  the  distance 
from  any  position  to  the  center  of  the  threat. 
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Figure  5.1:  Terrain  with  threats  formulated  as  an  exponential  function. 


Real  terrain  data  was  acquired  from  the  United  States  Geological  Survey  to 
incorporate  into  this  model[30].  The  data  was  found  in  tabular  format  relating  the  altitude 
to  the  locations  longitude  and  latitude,  with  data  points  spaced  approximately  every  48 
feet.  This  data  was  then  converted  to  matrix  form,  from  which  it  could  then  be  used 
as  / (x,  y ) .  Because  of  the  distance  between  the  sampled  altitude  points  in  the  matrix,  the 
data  was  then  smoothed  to  appear  more  continuous  and  to  remove  discontinuities  in 
altitude.  The  gradients  of  this  matrix,  along  both  the  x  and  y  directions,  were  calculated 
numerically  to  form  matrices  representing  fx  ( x ,  v)  and  /  (x,  y) .  The  gradients  of  these 

two  matrices  yielded  matrices  for  fxx(x,y)  ,  /vv(x,y)  and  /vv(x,y) . 
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Terrain  Elevation 


Figure  5.2:  Terrain  plot  of  an  area  near  Columbus  Ohio. 

For  this  portion  of  the  testing,  it  was  decided  to  use  a  section  of  terrain  near 
Columbus,  Ohio.  A  profile  of  this  terrain  can  be  seen  in  Figure  5.2.  In  this  graph,  the  x 
and  y-axes  depict  the  position  coordinates,  measured  in  feet,  such  that  the  x-axis  point 
north  and  the  y-axis  points  east.  The  altitude  of  the  terrain  is  measured  along  the  z-axis 
and  is  also  given  in  feet.  This  plot  depicts  a  square  plot  of  land,  with  10,000  feet  to  a 
side.  The  measurements  along  the  x  and  y-axes  are  relative  to  a  set  origin. 

5.2  Pseudo  3-D  and  3-D  results 

Optimal  trajectories  from  the  pseudo  3-D  formulations  that  include  wind  effects, 
moving  targets,  interior  point  constraints  in  the  form  of  waypoints,  and  moving  threats, 
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and  those  from  the  3-D  formulation  are  found  in  the  chapters  3  and  5  of  the  thesis  [19]. 
The  simulations  were  performed  with  generic  terrain  models  in  Figure  5.1  and  the  real 
terrain  data  in  Figure  5.2.  In  this  report,  we  compare  the  trajectories  obtained  from  both 
pseudo  3-D  and  3-D  reduced-order  formulations.  For  a  given  terrain,  a  flat  plane  with  a 
single  hill,  the  optimal  trajectories  found  to  navigate  it  for  both  minimum  time  and  terrain 
masking  flight  are  compared.  The  formulations  considered  include  the  simplified  and 
local  tangent  plane  equations  of  motion  for  the  pseudo-3D  case  as  well  as  the  constant 
velocity  and  varying  velocity  3D  equations  of  motion.  This  is  repeated  for  three 
different-steepness  hills. 

Figures  5.3  -  5.5  contain  the  results  for  the  minimum  time,  K  =  0,  case.  Each 
figure  portrays  the  results  from  a  different  hill  height.  The  top  two  plots  are  a  3D  view 
and  an  overhead  view  of  the  trajectories  from  the  pseudo-3D  case  while  the  bottom  two 
plots  depict  the  paths  for  the  3D  cases.  For  these  cases,  the  trajectories  from  using  the 
simplified  equations  of  motion  and  from  using  the  local  tangent  plane  equations  of 
motion  are  the  same  and  are  represented  by  the  black  line.  In  the  bottom  two  plots,  the 
black  line  represents  the  simplified  the  constant  velocity  3D  equations  of  motion. 


*  (M  0  0  y  (feet) 
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Figure  5.3:  Trajectories  for  K  =  0  and  hill  height  =  30. 


The  red  line  is  for  the  varying  velocity  3D  equations  of  motion.  The  plots  for  the  3D 
equations  of  motion  trajectories  are  the  same  as  depicted  in  Chapter  5  in  [19],  and  are 
repeated  here  for  convenience.  In  the  first  two  cases,  with  hill  heights  of  30  and  40  feet, 
the  constant  velocity  3D  trajectories  are  the  same  as  the  pseudo-3D  formulations; 
however,  in  the  steepest  hill,  the  3D  trajectory  begins  to  veer  around  the  hill.  In  all  three 
cases,  the  varying  velocity  trajectory  veers  around  the  hill  to  some  extent. 


X  (f8el)  °  0  y  (feet) 
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Figure  5.4:  Trajectories  for  K  =  0  and  hill  height  =  40. 


Some  details  from  the  various  cases  portrayed  in  these  figures  are  tabulated  below 
in  Table  5.1.  Here,  the  time  for  the  flight  and  the  cost  of  the  flight  can  be  compared  for 
each  set  of  equations  of  motion  investigated.  In  this  table,  it  can  be  seen  that  the  flight 
time  does  not  change  with  the  simplified  equations  of  motion,  while  the  time  increases 
for  the  local  tangent  plane  and  3D  equations  of  motion.  This  is  because  the  time  needed 
to  fly  vertically  is  ignored  in  the  simplified  equations  of  motion.  This  time  is  better 
accounted  for  with  the  local  tangent  plane  equations  of  motion,  but  these  final  times  are 
still  slightly  less  than  the  final  times  with  the  3D  equations  of  motion,  especially  when 
the  hill  steepness  is  greater.  The  final  times  for  the  varying  velocity  cases  are 
significantly  greater  due  to  the  loss  in  velocity  during  the  flights. 
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x  (feet) 


0  0 


y  (feet) 


*  (feel)  0  0  y  (feet) 


Figure  5.5:  Trajectories  for  K  =  0  and  hill  height  =  50. 


Hill  Height  =  30 

Hill  Height  =  40 

Hill  Height  =  50 

tf 

cost 

tf 

cost 

tf 

cost 

simplified 

7 

7 

7 

7 

7 

7 

local  tangent  plane 

7.06 

7.06 

7.1 

7.1 

7.15 

7.15 

3D 

7.06 

7.063 

7.11 

7.116 

7.16 

7.165 

varying  velocity 

7.27 

7.289 

7.37 

7.371 

7.45 

7.451 

Table  5.1  Trajectory  data  for  K  =  0  flights 
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Figures  5.6  -  5.8  contain  the  plots  for  the  results  for  if  =1.  These  plots  are  for 
cases  parallel  to  those  shown  in  Figures  5.3  -  5.5.  The  same  three  hill  heights  are 
displayed  with  the  four  sets  of  trajectories  depicted.  As  in  the  earlier  plots,  the  top  plots 


Figure  5.6:  Trajectories  for  K  =  1  and  hill  height  =  30. 

show  the  trajectories  from  using  the  simplified  and  local  tangent  plane  equations  of 
motion.  The  bottom  two  plots  depict  the  results  from  Chapter  5  in  [19]  for  the 
trajectories  using  the  3D  constant  velocity  and  varying  velocity  equations  of  motion.  In 
each  of  the  four  cases,  for  each  of  the  three  hill  steepness,  the  trajectories  appear  the 
same.  The  trajectory  always  curves  around  the  given  hill. 
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Table  5.2  contains  the  final  time  and  cost  information  for  the  trajectories  from 
each  set  of  equations  of  motion  and  for  each  hill  height  for  the  K  =  1  fonnulation.  It  can 
be  seen  here  that  there  is  very  little  difference  in  these  results  regardless  of  the  equations 
of  motion  used  or  the  hill  steepness. 


Figure  5.7:  Trajectories  for  K=  1  and  hill  height  =  40. 
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Figure  5.8:  Trajectories  for  K  =  1  and  hill  height  =  50. 


Hill  Height  =  30 

Hill  Height  =  40 

Hill  Height  =  50 

tf 

cost 

tf 

cost 

tf 

cost 

simplified 

8.42 

26.349 

8.42 

35.1314 

8.42 

43.9143 

local  tangent  plane 

8.42 

26.348 

8.42 

35.1305 

8.42 

43.9132 

3D 

8.42 

26.349 

8.42 

35.132 

8.42 

43.914 

varying  velocity 

8.41 

26.306 

8.42 

35.108 

8.42 

43.88 

Table  5.2  Trajectory  data  for  K  =  1  flights 


5.3  Six  DOF  optimal  trajectories  using  GESOP 
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While  reduced-order  optimal  trajectories  are  obtained  by  indirect  methods, 
GESOP  utilizes  a  direct  method.  Therefore,  it  is  first  verified  that  optimization  results  of 
GESOP  are  comparable  with  the  results  the  optimizations  in  [19]  by  implementing  the 
simplified  equations  in  Section  2.2  into  GESOP  with  the  same  boundary  conditions 
enforced.  It  is  shown  that  the  resulting  trajectories  are  almost  identical  (see  Chapter  5  in 
[18]).  This  ensures  that  numerical  methods  employed  for  optimization  do  not  generate 
differences.  The  first  optimizations  were  the  flights  around  one,  three  and  five 
mountains.  The  3DoF-Optimizations  were  used  for  the  optimizer  verification  presented  in 
Appendix  B  in  [18].  Later  when  the  6DoF  model  was  working,  the  same  optimizations 
were  done  with  it  to  investigate  its  behavior  for  flights  over  mathematical  formed  terrain 
as  in  (5.1). 

For  the  6DoF-Optimization  K  was  set  to  1  and  Terrain-Following  was  activated. 
Since  the  ground  under  the  optimal  path  is  almost  perfectly  flat  a  constant  absolute 
altitude  must  be  the  result.  As  shown  in  Figure  5.12  this  is  almost  perfectly  the  case 
except  a  very  small  variation  of  about  7ft  on  a  flight  path  length  of  more  than  1600ft. 
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Figure  5.9.  Trajectories  for  3  Mountains  -  3  DOF  GaTech  GESOP 
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In  figure  5.9  practically  identical  results  of  the  GaTech-Calculation  (dashed)  and 


the  GESOP-3DoF-Optimization  (dotted)  are  shown  with  the  corresponding  mountains. 


The  x-v- units  are  in  [feet]  here.  The  zoomed  trajectories  are  shown  in  the  right  figure  in 


Figure  5.9. 


As  a  next  step,  case  studies  with  different  operating  conditions  and  cost  functions 
are  extensively  carried  out  to  understand  the  influence  of  different  initial  conditions, 
dynamic  constraints,  and  different  cost  functions.  Table  5.3  illustrate  cases  studied  with 
the  full  6  DOF  GESOP  optimizations,  whose  definitions  for  case  details  can  be  found  in 
[18]. 


Case 

Terrain 
Masking 
optimization 
(K=1 ) 

Explicit 

Time 

optimization 

(K=0) 

Terrain 

following 

on 

Visibility 

Check 

High 

Penalty  for 
flying  under 
Terrainlevel 

Final  altitude 
constraint 

Comparable 

weighting 

Flight 

altitude 

Maximal 

Throttle 

1 .  Reference  1 

X 

X 

High  (50m) 

1.0 

2.  Reference  II 

X 

X 

X 

X 

X 

High  (50m) 

1.0 

3.  Reference  1,  Lower  flight 
altitude 

X 

X 

Low  (18m) 

1.0 

4.  Reference  1,  Lower  flight 
altitude  +  penalty  for  flying 
under  terrain  level 

X 

X 

X 

Low  (18m) 

1.0 

5.  Reference  1,  Half  Thrust 

X 

X 

High  (50m) 

0.5 

6.  Reference  1,  deactivated 
engine  moments 

X 

X 

High  (50m) 

1.0 

7. Reference  II,  Optimal 
Time 

X 

X 

X 

X 

X 

High  (50m) 

1.0 

8. Reference  II.  Visibility 
Check 

X 

X 

X 

X 

X 

High  (50m) 

1.0 

9. Reference  II.  Visibility 
Check  +  Terrain  following 

X 

X 

X 

X 

X 

X 

High  (50m) 

1.0 

10.  Reference  1  +  Fina 
Constraint 

X 

X 

X 

High  (50m) 

1.0 

11.  Reference  1,  Optima 
Time 

X 

X 

High  (50m) 

1.0 

12.  Reference  II,  Lower 
flight  altitude  +  penalty  for 
flying  under  terrain  level 

X 

X 

X 

X 

X 

Low  (18m) 

1.0 

13.  Flight  over  terrain  foi 
pop  up  threat  (threat  not 
active) 

X 

X 

X 

X 

X 

Low  (20m) 

1.0 

14.  Flight  over  terrain  foi 
pop  up  threat  (threat 
active) 

X 

X 

X 

X 

X 

Low  (20m) 

1.0 

Table  5.3  Parametric  Studies-Overview 
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5.4  Comparisons  between  Pseudo  3-D  results  and  GESOP  solutions  with 
Simulated  hills 

5.4.1  Flight  around  three  mountains 

Figure  5.10  shows  the  dotted  GESOP-pseudo  3-D  Optimization,  labeled  as 
“GaTech  3DOF-Opt”,  now  in  comparison  with  the  trajectory  followed  by  the  aircraft  in 
the  f-wing  simulator,  labeled  as  “GaTech  Sim”,  and  the  solid  6D0F-GESOP- 
Optimization,  labeled  as  “GESOP  6DOF-Opt”. 


2000 


-  -  GESOP  3DoF-Opt 
-  GaTech  Sim 


Figure  5.10  3  Mountains-3  DoF-GESOP  vs  Figure  5.10  3  Mountains-3  DoF-  with  and 
GaTech-Sim  &  GESOP  6  DoF  without  Turn  rate  constraint  and 

GESOP  6  DoF 


The  GaTech- Simulation  follows  the  commanded  trajectory  very  closely  up  to  the  sharp 
turn.  There  it  follows  the  command  with  a  short  delay  and  continues  on  a  path  parallel  to 
the  commanded  trajectory  with  a  constant  shift  until  the  end  of  the  simulation.  The  solid 
6DoF-Optimization  starts  and  ends  at  the  same  point  like  the  3DoF-Optimization. 
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This  is  not  surprising  because  these  points  are 
set  as  initial  and  final  boundary  conditions.  The 
trajectory's  course  is  also  about  the  same.  But 
the  dynamics  of  the  airplane  model  seem  to 
cause  a  smoothening  of  the  trajectory.  The  turn 
is  less  sharp,  it  begins  earlier,  ends  later  and  the 
maximal  extend  to  the  left  is  smaller.  It  looks 
like  the  maximal  turn  rate  is  restricted  by  the 
airplane's  dynamics.  So  an  interesting  question 
is  how  the  3DoF-Optimization  would  change  if  a  maximal  turn  rate  constraint  would  be 
set.  The  result  is  shown  in  figure  5.11.  The  dotted  line  is  the  3DoF-Optimization  as 
before,  the  solid  line  the  6DoF-Optimization  as  before,  and  the  dashed  line  the  new 
3DoF-Optimization  with  the  turn  rate  constraint.  The  3DoF-  and  6DoF-solutions  are 
almost  identical  now.  So  for  flights  over  flat  terrain  at  a  constant  altitude  a  turn  rate 
restriction  in  the  simple  3DoF-Simulation  can  produce  a  solution  very  similar  to  the 
6DoF-solution.  Most  probably  the  autopilot  of  the  simulator  would  have  minor  problems 
to  follow  the  commanded  trajectory,  too,  because  no  impossible  turn  forcing  the  airplane 
to  leave  the  commanded  trajectory  would  be  included  any  more. 


Figure  5.12:  3  Mountains  -  6DoF  - 
time-altitude 


5.4,2  Flight  around  Five  mountains 

As  can  be  seen  in  the  figures  5.13  and  5.14  the  3DoF-Optimizations  from  GaTech 
and  GESOP  are  identical  again,  in  figure  5.15  gets  clear  again  that  the  GaTech- 
Simulation  follows  with  some  delay  and  a  persisting  deviation  while  the  6DoF- 
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Optimization  follows  closely  with  a  slightly  lower  maximal  turn  rate.  In  figure  5.16  the 
maximal  turn  rate  of  the  dashed  3DoF-Gesop-Optimization  was  set  to  a  value  so  that  the 
3DoF-Result  is  similar  to  the  6DoF-Trajectory.  But  at  the  first  little  turn  it  is  obvious  that 


Figure  5.13:  5  Mountains  -  3DoF  GaTech-Gesop  Figure  5.14  1:  5  Mountains  -  3DoF  GaTech- 
-  xy  a  Gesop  -  xy  b 


Figure  5.15:  5  Mountains  -  3DoF  Gesop  vs.  Figure  5.16:  5  Mountains  -  3DoF  with/-out  turn 
GaTech  Sim  &  Gesop  6DoF  -  xy  constraint  vs.  6DoF  Gesop  -  xy 

for  a  good  turn  rate  constraint  approximation  it  must  be  not  just  constant,  but  velocity 
dependent.  The  limited  solution  has  a  bigger  distance  to  the  unlimited  solution  and  the 
6DoF-Optimization  is  able  to  follow  that  very  closely.  As  shown  in  figure  5.17  the  6DoF- 
velocity  is  still  low  there.  In  the  second  turn  the  6DoF-Simulation  takes  a  wider  way 
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around  the  mountain  and  sticks  closer  to  the  limited  solution  due  to  its  higher  velocity.  So 
a  velocity  dependent  turn  rate  restriction  in  the  3DoF-Optimizations  for  flights  in 
constant  absolute  altitude  would  approach  the  6DoF-Optimization  more  closely  than  the 
constant  restriction. 


time  [s] 


Figure  5.17:  5  Mountains  -  6D0F  -  velocity  u  over  time 

5.4,3  Flights  With  Moving  Targets  and  Threats 

In  addition  to  the  comparisons  of  the  3DoF-Target-Threat-Optimizations,  now  the 
6DoF-Optimizations  and  the  3DoF-Simulations  will  be  compared  to  the  3DoF- 
Calculations.  As  3DoF-Calculation  always  the  GaTech  solution  is  shown. 
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Figure  5.18  shows  the  6DoF-Optimization  in 
comparison  with  the  GaTech-Calculation. 
There  is  no  big  surprise,  the  6D0F- 
Optimization  takes  the  turn  a  little  smoother, 
but  this  is  the  only  difference.  In  the  GaTech- 
Simulation  case  in  figure  5.19  the  path  is 
followed  almost  perfectly  on  the  straight  line 
up  to  the  turn.  There  a  shift  is  generated  and 
doesn't  get  completely  eliminated  till  the  end 
point.  It  is  anticipated  that  an  I-factor  in  the 
autopilot  is  possibly  too  small  there. 

A  completely  different  course  than  the 
GaTech  calculation  takes  the  6DoF-Simulation 
in  the  1  threat  case  (figure  5.20).  Instead  of 
passing  in  front,  it  passes  behind.  But 
examinations  showed  that  this  pass-behind 
solution  can  also  be  produced  with  the  Gesop- 
3DoF-Optimization.  Simply  the  heading  for 
the  initial  guess  needed  to  be  changed  and  a 
solution  very  similar  to  the  6DoF-Optimization 
(dotted  line)  was  the  result.  But  it  wasn't 
possible  to  do  it  the  other  way  round  so  that  the 


I'XI  HD  JUC  *OQ  SX  SCO  M)  000  *D  tOO] 

•  1*1 


Figure  5.  18:  Moving-Target:  GaTech-Calc 
vs.  6DoF-Gesop  -  No  Threat 


wo 


1  MO  JX  XT  *E  5X  rtX  in  n  *0  HOC 

•1*1 

Figure  5.19.:  Moving-Target:  GaTech-Calc 
vs.  GaTech-Sim.  -  No  Threat 


■1*1 


Figure  5.20:  Moving-Target:  GaTech-Calc  vs. 
6DoF-Gesop,  3DoF-Gesop  restricted  -  One 
Threat 
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6DoF-case  would  follow  the  original  3DoF-case. 


Of  course  the  GaTech-Simulation's  autopilot  in  figure  5.21  follows  the 
commanded  trajectory  but  here  the  gap  between  the  commanded  and  the  simulated  path 
is  wider. 


Figure  5.21:  Moving-Target:  GaTech-Calc  Figure  5.22:  Moving-Target:  GaTech-Calc 
vs.  GaTech-Sim.  -  One  Threat  vs.  6DoF-Gesop,  3DoF-Gesop  restricted  - 

Two  Threats 

It  closes  faster  than  before,  but  not 
completely  till  the  end  as  well.  The  second 
threat  doesn't  change  the  situation  in 
general  any  more.  Also  here  the  6DoF- 
Optimization  passes  behind  both 
mountains,  the  original  3DoF- 
Optimizations  pass  in  front,  and  again  a 
3DoF-Optimization  with  a  different 


wo  an  an  «o  soo 


TO  SC  TO 


Figure  5.23:  Moving-Target:  GaTech-Calc 
vs.  GaTech-Sim.  -  Two  Threats 


heading  for  the  initial  guess  produces  a  solution  similar  to  the  6DoF-version.  The 
autopilot  in  figure  5.23  reacts  as  described  before,  but  it  should  be  noticed  that  the  gap 
after  the  second  turn  is  bigger  than  after  the  first  one.  This  could  be  because  the  second 
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turn  is  simply  harder,  but  at  least  a  portion  is  owed  to  the  fact  that  the  gap  isn't  closed  at 
the  beginning  of  the  second  turn.  Therefore  more  turns  could  lead  to  problems  with  the 
stability  of  the  simulation  when  it  gets  too  far  off  the  commanded  track. 


5.5  Comparisons  between  Pseudo  3-D  results  and  Gesop  solutions  with 
the  real  terrain 


5.5.1  Terrain-Masking  (K=l) 

The  first  3DoF-6DoF  comparisons 


over  real  terrain  are  optimizations  with 
K=1  (Terrain-Masking).  The  Gesop- 
Optimization  is  the  lower  case  of  the 
comparison  with  K=1  Low-Flight  -  High- 
Flight  shown  above.  In  figure  5.24  the 
GaTech  (black)  -Calculation  (dotted),  its 
simulation  (solid),  together  with  the 
Gesop-Optimization  (white)  are  drawn.  It 


Figure  5.24:  Ref.II,  K=1  vs.  GaTech-Calc  & 

can  be  seen  immediately  that  the  GaTech  -Sim.-3D 

solution  takes  a  completely  different  way  compared  to  the  Gesop  trajectory.  Despite  the 
fact  that  the  optimizer  with  K=1  is  set  to  fly  over  terrain  which  is  as  low  as  possible,  the 
GaTech-Calculation  takes  a  path  over  the  yellow  ridge.  In  addition  the  simulation  has 
significant  problems  to  follow  the  commanded  path  (figure  5.25).  This  even  leads  to  a 
path  where  the  airplane  would  crash  into  the  terrain  at  the  first  hill  at  this  altitude  (break 
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in  the  solid  black  line  in  figure  5.25).  Also  the  commanded  final  altitude  cannot  be 
reached.  The  reason  for  that  is  most  probably  the  fact  that  the  simulator  tries  to  fly  a 
constant  velocity.  A  steeper  decent  would  accelerate  the  airplane  over  the  commanded 
constant  velocity,  so  the  Autopilot  descends  slower.  Regarding  the  flight  path  length  the 
Gesop-Optimization's  path  is  a  little  shorter  (figure  5.26).  Also  the  integral  of  the  terrain 
altitude  over  the  path  length  shows  the  significant  quality  difference  between  the  paths. 


Figure  5.25:  Refill,  K=1  vs.  GaTech-Calc  &  -  Figure  5.26:  Ref.II,  K=1  vs.  GaTech-Calc  &  - 
Sim.  -  path-altitude  Sim.-xy 

The  GaTech  path  integral  with  1 15,838.64  is  about  34.2%  higher  than  the  Gesop  path 
integral  with  86,308.24. 


5.5.2  Optimization  of  time  (K=0) 

Now  instead  of  using  K=1  for  Terrain-Masking  optimization,  K=0  for  Time- 
Optimization  was  used  in  figures  5.27-5.29.  The  Gesop  case  this  time  is  the  higher  case 
of  the  Time-Optimization  shown  above. 
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Compared  to  the  previous  case,  the 
GaTech-Optimization  again  takes  a 
different  path  than  the  6DoF-Gesop 
solution.  But  an  optimal  solution  with 
(unfortunately  not  comparable)  different 
weightings  on  altitude  change  cost  takes 
the  same  valley  and  sticks  much  closer  to 
the  GaTech  solution.  These  results  are 
shown  in  the  figures  5.30-5.32  below. 


Figure  5.28:  Ref.II,  K=0  vs.  GaTech-Calc  &  - 
Sim.-xy 


Figure  5.27:  Ref.II,  K=0  vs.  GaTech-Calc  & 
-Sim.-3D 


Figure  5.29:  Ref.II,  K=0  vs.  GaTech  -Calc  & 
Sim.  -  path-altitude 
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The  reason  that  despite  the  fact  that  for 
K=0  the  GaTech  Trajectory  is  not  a 
straight  line  is  the  use  of  the  pseudo  3-D 
formulation.  They  let  the  optimizer  realize 
that  flying  over  a  hill  lengthens  the  flight 
path  but  still  doesn't  tell  it  the  differences 
in  velocity  caused  by  climbs  /  descents. 


Figure  5.30:  Ref.I,  K=0  vs.  GaTech  -Calc  & 
-Sim.-3D  -  unequal  weighting 


Figure  5.31:  Ref.I,  K=0  vs.  GaTech-Calc  &  -  Figure  5.32:  Ref.II,  K=0  vs.  GaTech  Calc  & 
Sim.-xy  -  unequal  weighting  Sim.-  path-altitude  -  unequal  weighting 


Therefore  the  optimizer  tries  to  avoid  high  slopes  and  at  the  same  time  it  tries  to  find  the 
shortest  possible  connection,  which  is  probably  the  reason  why  it  flies  so  close  nearby  the 
red  top  mountain  at  a  relatively  constant  altitude  without  using  the  valley  to  accelerate 
and  reaching  the  destination  quicker.  Again  the  constant  velocity  constraint  inhibits  the 
simulator  to  be  able  to  reach  the  destination  at  the  commanded  altitude  and  it  is  not  able 
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to  follow  the  altitude  commands  closely  along  the  trajectory  (figure  5.32).  On  the  other 
side  the  x-v-command  is  followed  very  well  like  before  (figure  5.31).  This  time  the 
simulation's  trajectory  does  not  cross  the  terrain,  but  if  the  starting  altitude  would  be 
lower  it  would. 

In  comparison  the  Gesop  solution  consequently  tries  to  fly  over  low  profile  as 
much  as  possible.  First  the  obligatory  decent  to  gain  speed  and  then  it  climbs  the  first  hill 
at  the  latest  position  possible  for  this  trajectory.  Afterwards  it  flies  directly  in  the  middle 
of  the  first  valley  to  descent  as  soon  as  possible  again.  So  the  average  overflown  altitude 
is  lower  in  Gesop  and  therefore  the  Flight  Path  -  Terrain  altitude  integral  with  86.308.24 
is  about  31,59%  smaller  than  1 13,574,47  for  the  GaTech-Optimization. 

5.5.3  Flight  over  terrain  for  a  Threat 

The  white  line  in  figure  5.33 
represents  the  6DoF-solution,  the  black 
dashed  line  the  GaTech-Optimization,  the 
black  solid  line  the  GaTech- Simulation 
when  a  threat  pops-up. 

In  figure  5.34  the  trajectories  in  the  x-y- 
plane  are  shown  and  in  figure  5.35  the 
altitude  over  the  flight  path  length  can  be 
seen.  Despite  the  fact  that  both 
optimizations  were  made  with  K=1  and  the  terrain  level  input  for  the  Lagrange  Cost 
function  was  normalized  to  the  same  values,  there  still  remains  a  significant  difference. 


/ 

I 


Figure  5.33:  Pop-Up-Threat  Gesop  (white) 
vs.  GaTech  (black)  - 
Calc.(dashed)  &  -Sim.(solid)  -  3D 
-  No  Threat 
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Figure  5.34:  Pop-Up-Threat  Gesop  vs.  GaTech-  Figure  5.35:  Pop-Up-Threat  Gesop  vs.  GaTech- 
Calc.  &  -Sim.  -  xy  -  No  Threat  Calc.  &  Sim.  -  path-altitude  -  No  Threat 

The  6D0F  solution  takes  a  longer  way  compared  to  the  GaTech  solution.  The  reason  for 
this  difference  is  most  probably  caused  by  the  advantage  the  6DoF-Model  draws  out  of 
the  fact  that  it  is  able  to  gain  speed  much  quicker  by  flying  through  the  blue  sink  and 
even  more  avoiding  the  higher  but  shorter  path  the  GaTech  solution  takes.  As  described 
earlier  a  constant  velocity  is  commanded  for  the  3DoF-simulation,  which  does  not  take 
into  consideration  the  extension  of  the  path  and  slow  down  by  first  going  up  and 
afterwards  going  down  a  hill.  For  a  flight  going  from  a  point  A  to  a  point  B  it  is  faster  to 
fly  through  an  imaginary  valley  rather  than  to  fly  over  an  imaginary  mountain  with  the 
same  size  because  the  average  speed  is  much  higher  when  accelerating  in  the  beginning 
of  the  maneuver.  The  6DoF-optimizer  realizes  the  advantage  of  a  shorter  flight  time  and 
therefore  smaller  Lagrange  Integral  by  quickly  gaining  speed  through  the  valley  in  the 
beginning  and  also  tries  to  avoid  higher  mountains  since  they  slow  down  the  airplane.  A 
look  at  figure  2.53  shows  that. 
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The  6D0F  flight  is  slightly  longer  but  the  average  flight  altitude  of  the  Gesop 
solution  is  much  lower  than  the  GaTech-Calculation  /  -Simulation.  The  altitude 
difference  at  the  end  is  most  probably  caused  by  the  different  interpolation  methods  used 
(Cubic-Spline -Int.  in  Gesop  vs.  Linear-Int.  for  GaTech).  So  at  this  example  the 
differences  between  3DoF  and  6DoF-Optimiations  come  out  very  clearly. 


5.5.4  Flights  with  a  Pop-UP  Threat 

Now  during  the  flight  presented  in  the  previous  paragraph  a  Pop-Up-Threat 
suddenly  appears  after  some  seconds  of  the  flight  on  the  optimal  trajectory  of  the  airplane 
so  that  the  airplane  must  find  a  new  trajectory  to  avoid  the  threat.  The  threat  itself  is 
again  represented  by  a  Gaussian  hill.  Figure  5.36  shows  the  situation.  The  solid  line 
shows  the  normal  flight  path,  at  the  point  where  the  dashed  line  comes  out  of  the  solid 
line,  the  wired  hill  representing  the  threat  suddenly  appears.  At  this  point  a  new 
optimization  is  started  with  all  actual  state -values  as  initial  boundary  conditions  for  the 
new  optimization.  The  optimizer  calculates  the  new  trajectory  which  is  represented  by 
the  dashed  line  in  Gesop.  The  GaTech-Solver  doesn't  need  to  use  these  initial  boundary 
conditions  since  it  just  creates  heading  commands  for  the  simulator  which  the  Autopilot 
tries  to  follow  then.  So  sudden  turns  are  allowed  there.  The  effect  can  be  seen  in  figure 
5.37.  The  dashed-dotted  line  is  the  heading  command,  the  solid  line  is  the  simulator's 
trajectory  trying  to  follow  the  heading  command.  Coming  from  the  right  bottom,  after  the 
appearance  of  the  Pop-Up-Threat  the  heading  command  suddenly  jumps  clockwise.  The 
airplane  of  course  can't  follow  such  quickly  and  tries  to  adapt  with  the  delay  it  needs  to 
turn. 
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Figure  5.36  Pop-Up-Threat-Gesop  Solution-3  Figure  5.  37  Pop-Up-Threat-GaTech- 
D-with  Threat  (white  solid)-without  Threat  Calc.(dashed)  &  -Sim  (solid)-with  Threat 
(white  dashed) 


The  6DoF-Optimization  creates  this  smooth  turn  already  during  the  optimization  because 
the  aircraft  model  with  its  dynamics  is  already  built  into  the  optimization  loop.  Since  the 
basic  trajectories  in  the  paragraph  before  were  too  different  to  produce  a  comparable 
result  with  the  same  (threat-)hill,  the  hills  for  the  two  optimizations  were  located  at 
different  positions  so  that  their  tops  intersect  with  the  corresponding  trajectory.  Then  the 
optimizations  with  the  Pop-Up-Threat  were  set  in  such  a  way  that  the  distances  of  the 
airplane  to  the  appearing  threat  were  about  the  same  in  both  cases.  In  figure  5.38  a 
comparison  between  the  GaTech-Simulation  (black)  and  the  Gesop-Optimization  (white) 
for  a  Pop-Up-Threat  appearing  very  close  in  front  of  the  airplane  is  shown.  It  is  the  same 
case  like  in  figure  5.36  and  5.37.  The  hills  were  left  out  to  keep  the  figure 
understandable,  but  figure  5.36  shows  the  same  case  with  the  hill  for  the  Gesop-part.  The 
points  where  the  threat  appears  are  marked  by  additional  small  circles  in  the  figures  5.39 
and  5.40.  The  behavior  of  both  solutions  is  very  similar.  Both  try  to  avoid  the  hill  by 
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flying  around  it.  The  fact  that  they  take  different  ways  around  the  hills  is  not  significant 
because  slight  hill  position  changes  can  already  switch  that. 

As  can  be  seen  in  figure  5.39  the 
turn  rates  in  the  GaTech-Simulation  seem 
to  be  a  little  higher  than  in  the  6D0F- 
Optimization.  Reasons  for  that  might  be 
the  penalty  on  the  control  derivatives,  or 
the  relatively  wide -meshed  control  grid, 
but  one  certain  reason  is  that  the  average 
velocity  in  the  GaTech-Simulation  was 
lOOft/s  (30.48m/s)  while  the  average  speed 
of  the  6DoF-Optimization  with  109.25  ft/s 
(33.3m/s)  was  about  10%  higher. 

In  contrast  to  the  undisturbed  flight  both  solution  also  fly  a  lower  profile  in  the 
case  with  threat  (figure  5.40).  The  explanation  for  that  are  most  probably  the  sharp  turns. 
In  the  6DoF-case  the  optimizer  tries  to  keep  the  velocity  as  high  as  possible  which  leads 
to  a  descent  during  the  sharp  turn  because  the  lift  decreases  and  a  higher  lift  during  the 
sharp  turn  would  result  in  a  higher  drag. 

On  the  other  hand  for  the  GaTech  case  where  the  simulator  tries  to  fly  a  constant 
velocity  the  loss  in  lift  by  the  sharp  turn  maybe  can't  be  compensated. 


Figure  5.38:  Pop-Up-Threat  appearance 
Gesop  (white)  vs.  GaTech 
Sim.(black)  Dashed:  With 
appeared  threat,  Solid:  No 
appeared  threat  -3D 
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Figure  5.39:  Pop-Up-Threat  appearance 
Gesop  (green)  vs.  GaTech  Sim.(black)  Dashed: 
With  appeared  threat,  Solid:  No  appeared 
threat -xy 


Figure  5.  40:  Pop-Up-Threat  appearance 

Gesop  (green)  vs.  GaTech  Sim.(black)  Dashed: 
With  appeared  threat,  Solid:  No  appeared 
threat  -  path-alt. 


5.5.5  Earlier  pop-up  Threat 

In  the  dotted  lines  in  figure  5.41  shows  the  result  when  the  Pop-Up-Threat 
appears  earlier  and  the  aircraft  has  more  time  to  avoid  it.  As  expected  the  trajectory 
leaves  the  original  path  earlier.  After  some  time  it  matches  the  dashed  trajectory  of  the 
first  Pop-Up  Threat  in  both  cases.  This  shows  that  behind  the  threat  the  optimal  path  is 
about  constant  and  can  be  reproduced  even  under  different  circumstances. 
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Figure  5.41:  Early  Popup-Threat  Gesop  (white)  vs.  GaTech  (black)  Calc.  &  Sim.  -  Dotted: 
early  appearance  of  threat,  Dashed:  late  appearance  of  threat  -  3D 
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6.  Flight  Demonstrations 


This  section  includes  flight  demonstration  results  in  which  an  optimal  trajectory  is 
implemented  in  real-time  using  the  reduced-order  formulation.  The  pop-up  obstacle  is 
simulated  as  a  pop-up  mountain  in  a  virtual  environment. 

6. 1  Integration  Details 

The  planner  was  integrated  into  the  Rmax  environment  for  execution  on  the 
second  flight  computer  (onboard2).  A  list  of  tasks  being  performed  on  onboard2  is  shown 
in  Figure  6.3. 


The  general  sequence  of  events  when  flying  the  planner  includes  the  following. 


1.  On  gcs,  load  obstacle  database  into  memory  on  the 
ground  station  (gcs) . 

a.  wdbset  gcswdb 

b.  wdbclear 

c.  wdbdatumned 

d.  @ropdb 

e.  sceneO . redraw=l 

f.  scenel . redraw=l 

g.  scene2 . redraw=l 

2.  On  gcs,  instruct  onboardl  to  start  sending  extmanO 
messages 

a.  rc  obDatalinkSet . sendExtManO  =  1 

3.  On  onboard2;  Set  planner  parameters;  an  important 
parameter  is,  hc_const,  which  must  be  set  to  the 
nominal  altitude  above  the  terrain  that  one  wants  the 
vehicle  to  fly: 

a.  cd  reducedOrderPlanner 

b .  W  =-0.05 
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c.  V_stall  =  5 

d.  flag_real  =  1 

e .  e  =1 

f.  w_end  =  500 

g.  hc_const  =  300 

h.  scale_lam  =  1 

i.  T  =  50 


4.  On  onbard2,  set  up  history  recording  by  executing 

a.  @rophist 

5.  On  onbard2,  get  ready  for  loading  the  obstacle  database 

a.  wdbset  ob2wdb 

b.  wdbclear 

6.  On  gcs,  issue  a  command  that  transmits  the  datum  used 
to  initialize  onboardl,  to  onboard2 . 

a.  wdbdatumob2 

7.  On  onbard2,  load  the  obstacle  database 

a.  @ropdb 

8.  On  onboard2.  Start  the  rop  planner  task 

a.  taskStart  rop 

9.  On  onboard2,  instruct  the  onboard2  main  loop  to  call 
ropTraj  when  an  external  maneuver  message  (extmanO)  is 
received 

a.  onboard2 . runRop=l 

The  above  steps  are  encoded  as  input  fdes  and  executed  in  the  following  sequence. 

1.  on  gcs  do  Sropgcs 

2.  on  ob2  do  @ropob2 

3.  on  gcs  do  wdbdatumob2 

4.  on  ob2  do  @ropob2b 


At  this  point,  the  planner  task  is  executing  on  onboard2.  Various  waypoints  may 
be  setup  on  the  gcs,  some  that  may  include  passing  through  obstacles.  The  various 
pseudo-codes  for  the  tasks  executing  on  onboard2  are  shown  in  Figure  6.3.  An  overview 
of  the  messages  being  send  and  received  across  the  various  computers  are  shown  in 
Figure  6.4. 

6.2  Simulation 
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As  an  example,  the  planner  was  executed  in  simulation.  In  Figure  6  1 :  Example  gcs 
view  during  before  using  the  planner  segment  1  passes  through  an  obstacle.  In  this  case  the 
obstacles  are  hills  with  a  maximum  altitude  of  300  ft. 


segment  1 

\ 

segment  2 

H  ■  Scene  Window  1 

■  show  plan  map  J 

■  fixed 

h| 

IKJ 

I  1  I  r~ BS 1  ■  1  ;'"j 

segment  0 


obstacles 


Figure  6  1:  Example  gcs  view  during  before  using  the  planner 


The  planner  running  on  onbard2  may  be  instructed  to  replan  segment  1  (between 
waypoint  1  and  waypoint  2),  by  issuing 


replan  1 


at  the  gcs  console.  Figure  6  2  now  shows  the  trajectory  for  segment  1  generated  by  the 
planner  in  blue. 
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re-planned  segment  1 

\ 

r 

■  Scene  Window  1 

□11]® 

show  plan...  |  map  | 

fixed 

Figure  6  2:  gcs  after  segment  1  has  been  replanned.  The  blue  segment  indicates  trajectory 
generated  by  the  planner 

Making  sure  that  waypoint  2  has  the  type  MAN_EXT,  the  trajectory  may  be 
uploaded  and  executed  using  the  following  commands. 


tra j  Upload 
tra j  Go 


Figure  6.5  shows  the  aircraft  at  various  instants,  tracking  the  planner  generated 
trajectory.  Figure  6.6,  shows  a  close  up  of  the  position  response  around  the  hill  and 
Figure  6.7  shows  the  position  response  for  the  entire  set  of  waypoints. 
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6.3  Experimental  Results 

A  flight  test  was  conducted  on  the  11th  of  July  2007  -  utilizing  the  GTYak  research 
UAV. 

Pictures  (files  starting  with  GTYak): 

http://uav.ae.gatech.edu/pics/f070711  flight  FormationFlight  ObstacleAvoidance/ 

Video: 

http://uav.ae.gatech.edu/videos/v07071  lb  1  planner  l.wmv 

During  the  test,  a  simulated  hill  was  placed  at  0,0,0  with  a  maximum  altitude  of  300  ft. 
The  setup  for  simulation  and  flight  is  exactly  the  same,  except  that  gcs,  onboard  1  and 
onboard2  are  executing  on  different  computers.  Figure  6.8,  shows  the  trajectory  around  a 
simulated  hill.  The  offset  in  lateral  position  is  due  to  the  initial  condition  offset  when 
entering  the  external  planner’s  segment.  This  offset  may  be  alleviated  by  introducing  a 
straight  line  segment  that  passes  through  the  last  waypoint  before  the  planner’s  segment 
is  executed.  Figure  6.9  shows  position  plots  versus  time. 

6.4  Relevant  Files 


File  name 

Description 

rop.h,  rop.cpp 

Standalone  planner  algorithm 

rop.db 

declares  data  structures  used  by  the  planner 
during  execution.  It  also  contains  data 
structures  used  to  store  planned  segments 

wdb.h,  wdb.cpp 

World  Database.  Contains  code  to 
represent,  add,  and  remove  obstacles.  Also 
contains  a  function  to  query  local  terrain 
altitude  taking  obstacles  into  account. 

wdb_scene.h,  wdb_scene.c 

World  Database.  Contains  code  to  draw 
obstacles  in  the  scene  window 

ropint.h,  ropint.cpp 

main  code  that  integrates  the  planner  into 
the  GTMAX  environment.  Contains 
procedures  to  generate  grid  data  from  the 
world  database,  makes  coordinate 
conversions  between  planner  and  GTMAX 
NED  coordinates.  It  also  contains  functions 
which,  given  a  pre-planned  trajectory 
segment,  will  use  interpolation  to  generate 
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trajectories  real-time  during  flight. 

ropgcs.inp 

execute  during  flight  test  at  gcs,  loads 
database  from  ropdb.inp,  sends  extmanO 

ropob2.inp 

some  parameters  for  the  planner,  activates 
ob2wdb,  the  onboard2  world  database  and 
clears  it,  also  setups  up  history  recording 
variables 

ropob2b.inp 

execute  on  ob2,  loads  database  from 
ropdb.inp  which  should  be  uploaded  to 
ob2,  also  starts  the  planner  thread 

ropdb.inp 

contains  list  of  obstacles  and  loaded  at  gcs 
and  ob2 

ropsim.inp 

setup  up  sim  environment  for  testing  rop 

rophist.inp 

history  variables  for  data  recording 
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onboard2 


Main  Thread  (updateOnboard2):  Planner  Thread  (ropTask): 


-  while  not  quit: 

Planner  Thread  (ropTask) :  ^ 

-  read  onboard2PortIPC 

-  while  not  quit: 

-  if  extmanO  was  received 

-  loop  through  each  segment: 

-  call  ropTraj 

-  if  segment  needs  replanning 

-  endif 

-  call  gengrid 

-  send  extmanl  (in  onboard2 . cpp) 

-  execute  planner 

-  endwhile 

-  if  valid  trajectory  generated 
-  set  validsolution  flag 

ropTraj : 

-  save  trajectory 

-  if  man  type  is  MAN  EXT  and  the  relevant 

-correct  coordinates 

segment's  validsolution  flag  is  true 

-smooth  trajectory 

if  first  time  through 

-  send  trajectory  to  gcs 

-  save  start  time  in  seg.t 

and  save  it  to  file 

if  seg.t  not  greater  than  seg.trajtf 

-  else 

-  use  getTraj  to  get  interpolated 

-  unset  validsolution  flag 

solution  at  trajectory  time  t 

-  endif 

-  fill  in  extmanl  with  trajectory 

-  endif 

information 

-  endloop 

-  else 

-  endwhile 

-  increment  extmanl .manlndex  to  move 

-  gengrid: 

to  next  waypoint 

-  figure  out  extents  of  database 

-  endif 

and  set  coordinate  corrections 

-  endif 

-  probe  world  database  to  get 
altitude  map 

-  use  central  differences  to 
generate  directional  derivatives 

-  save  grid  to  file 

Figure  6.3:  Overview  of  tasks  on  onboard2 
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onboardPortDatalinkl 

onboardPortDatalink2 


onboard  1 

onboardPortIPC 

onboard2PortIPC 

onboard2 

rx-gcs : obDatalinkMessageRopReplan 
tx-ob2 : obDatalinkMessageRopReplan 
rx-ob2 : obDatalinkMessageExtTra j [0-19] 
tx-gcs : obDatalinkMessageExtTra j  [0-19] 
tx-ob2 : obDatalinkMessageExtManO 
rx-ob2 : obDatalinkMessageExtManl 


rx-obl : obDatalinkMessageRopReplan 
tx-obl : obDatalinkMessageExtTra j [0-19] 
rx-obl : obDatalinkMessageExtManO 
tx-obl : obDatalinkMessageExtManl 


tx-obl : gcsDatalinkMes sageRopReplan 
rx-obl : gcsDatalinkMessageExtTraj [0-19] 

Figure  6.4:  Overview  of  transmission  and  reception  of  messages  relevant  to  the  planner 
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Figure  6.5:  Images  at  various  instants  of  the  maneuver 
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Figure  6.6:  Plot  of  trajectory  around  hill.  extmanO  represents  the  actual  trajectory  of  the  vehicle, 
while  extmanl  represents  the  commanded  trajectory 
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Figure  6.7:  Plot  of  the  various  position  command  (extmanl)  and  response  (extmanO)  versus  time. 
The  segment  where  man  type  is  6  corresponds  to  the  segment  generated  by  the  planner. 
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Figure  6.8:  Plot  of  trajectory  around  a  simulated  hill. 


80 


time  Is 


time  Is 


time  Is 


time  Is 


Figure  6.9:  Plot  of  position  responses  during  flight,  extmanl  represents  the  command  and  extmanO 
represents  the  response.  The  segment  where  mantype  is  6  corresonds  to  the  trajectoyr  generated 
by  the  planner 
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7.  Conclusions  and  Recommendations 


In  the  Phase  II  effort  we  extensively  addressed  trajectory  optimizations  in  incrementally 
complex  formulations  starting  from  the  pseudo  3-D  formulation,  via  3-D  point  mass  model  to  the 
full  6  DOF  aircraft  dynamics.  Optimizations  over  mathematical  and  real  terrains  were 
investigated  with  moving  targets,  moving  threats,  and  pop-up  threats  while  the  aircraft  moves 
along  the  terrain.  The  pseudo  3-D  reduced-order  formulation  is  flight  demonstrated  for 
replanning  an  optimal  path  when  a  virtual  obstacle  pops  up  using  a  fixed  wing  flight  test-bed. 
This  illustrates  that  the  resulting  optimal  trajectory  from  the  pseudo  3-D  fonnulation  is 
numerically  efficient  and  can  be  implemented  for  trajectory  optimization  for  autonomous  aerial 
vehicles. 

Comparison  studies  of  the  pseudo  3-D  solutions  with  the  optimal  solutions  obtained 
using  the  full  6  DOF  aircraft  dynamics  reveal  various  subtle  issues  in  path  planning  regarding 
the  sub -optimality  of  the  reduced-order  solution.  In  the  case  of  terrain-masking  and  flight  time 
minimization,  our  research  during  the  Phase  II  efforts  led  to  the  following  observations. 

First,  towards  the  task  of  terrain-following,  the  pseudo  3-D  fonnulation  is  associated  with 
the  reduced-order  point-mass  equation  that  incorporates  the  terrain  information  in  the  phase  of 
problem  formulation.  This  leads  to  an  trajectory  that  maintains  a  constant  altitude  over  the 
terrain  map  with  a  constant  maximal  velocity  to  minimize  the  flight  time[13].  While  this  is  a 
reasonable  assumption  for  those  trajectories  in  which  vertical  portion  of  velocities  are  negligible, 
this  assumption  may  be  violated  when  the  aircraft  flies  over  a  series  of  steep  hills  and  deep 
valleys.  In  case  of  6  DOF  optimizations,  terrain  following  is  achieved  by  adding  additional  cost 
tenns  in  the  perfonnance  index.  Figure  7.1  compares  each  costs  in  the  case  of  terrain-following 
(K=l)  for  the  optimal  trajectory  with  the  real  terrain  data,  which  reveals  that  the  cost  for  the 
terrain-following  dominates  the  overall  performance  index.  Figure  7.2  shows  the  same  costs  as  in 
Figure  7.2  for  a  non-optimal  solution  with  the  6  DOF  aircraft  dynamics.  Figure  7.2  reveals  that 
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there  is  a  flight  trajectory  for  which  additional  costs  for  maintaining  a  constant  altitude  over  the 
terrain  is  at  the  same  magnitude  of  the  terrain  following  (at  least  for  a  while),  and  this  implies 
that  the  reduced-order  optimal  solution  that  minimizes  the  only  portion  of  the  terrain-following 
cost  may  quite  deviate  from  the  6  DOF  trajectory. 


Figure  7.1:  Comparison  of  costs  for  optimal  Figure  7.2:  Comparison  of  costs  for  a  non¬ 
trajectory  for  terrain- folio  wing  optimal  trajectory  for  terrain-  following 

The  spikes  shown  in  Figure  7.2  also  imply  that  the  rate  of  control  signals  should  be  constrained 
to  result  in  a  flyable  trajectory.  However,  this  has  to  induce  additional  constraints  on  the  states 
used  as  controls  in  the  reduced-order  formulation.  This  is  also  related  to  the  next  issue. 

Second,  heading  constrains  with  the  reduced-order  formulation  leads  to  very  close 
optimal  trajectory  to  the  6  DOF  optimal  solution  in  case  of  generic  hills  in  which  optimal 
trajectories  veer  around  hills  (for  example  see  Figure  5.10  and  5.16).  With  the  real  terrain  data, 
the  main  differences  between  the  pseudo  3-D  solutions  and  6  DOF  optimal  solutions  result  from 
altitude  variations.  As  a  matter  of  fact,  the  pseudo  3-D  solutions  are  not  flyable  in  the  f-wing 
simulator  in  some  cases.  Considering  that  heading-constraint  generates  a  trajectory  close  to  an 
optimal  trajectory  with  6  DOF  optimizer,  this  observation  implies  that  additional  constraint  or 
formulation  is  required  to  address  neglected  part  of  the  vertical  dynamics. 
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Third,  it  was  found  that  an  optimal  solution  may  be  sensitive  to  initial  conditions.  For 
example,  in  Figure  5.20  in  Section  5.4.3,  the  pseudo  3-D  solution  passes  a  threat  in  front  while 
the  6  DOF  optimal  solution  passes  the  target  behind.  Changing  an  initial  guess  for  the  heading 
angle  resulted  in  a  solution  that  is  close  to  the  6  DOF  optimal  trajectory.  An  interesting  question 
is  how  different  costs  are  for  these  seemingly  discrete  trajectories.  The  study  reveals  that  the 
difference  in  cost  is  not  significant.  Likewise,  for  6  DOF  optimization,  several  different 
trajectories  of  the  same  level  of  costs  were  available,  and  small  changes  in  the  weighting  factors 
in  the  performance  index  could  lead  to  quite  different  trajectories  as  a  result  of  more  possible 
combinations  for  states  and  control  variables.  Following  the  above  observations,  further  research 
for  on-line  path  optimization  is  recommended  to  address  the  following. 

First,  in  the  current  effort,  the  pseudo  3-D  optimization  led  to  similar  trajectories,  with 
comparable  cost,  to  those  of  3-D  formulation  as  shown  in  Section  5.2.  This  became  a  basis  for 
ensuing  study  of  extending  the  reduced-order  formulation  into  multi-vehicle  fonnulation  that 
will  benefit  the  research  on  multi-vehicle  fonnation  and  cooperation  in  performing  tasks.  More 
quantitative  analysis  is,  however,  required  with  varying  degree  of  complex  terrains.  In  particular, 
analysis  of  the  degree  of  the  sub-optimality  due  to  neglected  portion  of  the  dynamics  will 
provide  clear  insight  on  the  feasibility  of  the  current  formulation  in  more  general  terrain  and 
mission  scenarios. 

Second,  a  main  drive  for  the  reduced-order  formulation  was  numerical  efficiency  in  a 
general  setting  of  solving  necessary  conditions  for  optimality  with  the  singular  perturbation 
theory  as  an  analysis  tool  for  the  neglected  portion  of  the  dynamics.  However,  formulating  the 
aircraft  dynamics  in  a  singular  perturbation  form  requires  identification  of  a  parasitic  parameter 
which  reveals  inherent  time  scales  in  the  aircraft  dynamics,  which  is  probably  most  difficult  part 
of  applying  the  singular  perturbation  method.  As  a  matter  of  fact,  a  large  effort  has  been  made 
for  a  systematic  methodology  for  selecting  time-scaled  state  variables  in  general  nonlinear 
optimal  control[31-34],  and  writing  the  original  dynamics  in  a  singular  perturbed  form  is  still  a 
main  theoretical  issue[14].  Perfonning  more  rigorous  study  on  the  fonn  of  the  terrain-following 
dynamics  may  allow  for  proper  corrections  if  the  reduced-order  solution  leads  to  unflyable 
trajectories  as  has  been  done  in  [17].  An  alternative  method  based  on  pseudo-control 
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hedging[35]  allowed  the  vehicle  to  fly  a  modified  trajectory  in  simulation,  which  is  discussed  in 
the  last  paragraph  of  this  section. 

Third,  finding  equivalent  constraints  in  the  reduced-order  formulation  in  case  when  the 
real  control  signals  are  position  and  rate  limited  is  worth  study.  In  particular,  for  the  current 
scenario  of  terrain-following,  the  assumption  of  constant  velocity  over  the  constant  altitude  over 
the  terrain  is  unrealistic  with  the  real  terrain  and  led  to  a  solution  that  differs  from  the  true  6 
DOF  optimization.  However,  properly  constrained  heading  rate,  the  resulting  trajectory  is  very 
close  to  the  full-order  optimizer.  This  implies  that  further  study  on  adding  constraints  on  the 
vertical  motion  may  lead  to  similar  results  as  was  done  for  X-Y  trajectories.  Some  directions, 
such  as  limiting  the  angle  of  attack  or  formulating  the  reduced-order  problem  in  modified  energy 
form,  were  proposed  in  Phase  I  final  report. 

Fourth,  while  the  pseudo  3-D  formulation  employed  in  this  study  was  an  attempt  to 
reduce  the  number  of  optimization  variables  by  exploiting  analyticity  of  optimally  conditions  in 
a  singular  perturbed  dynamics,  various  schemes  that  tries  to  reduce  the  number  of  optimization 
variables  also  exist.  For  real-time  trajectory  generation,  the  notion  of  differential  flatness  has 
recently  been  introduced  and  attracted  many  researchers  because  the  reduction  of  the  number  of 
variables  can  be  performed  in  utilizing  inherent  structure  of  the  dynamical  systems  in  flat 
systems,  and  therefore  dynamic  constraints,  the  equation  of  states,  are  automatically  met[36-38]. 
While  this  sounds  very  promising  in  tenns  of  numerical  computation,  recent  research  reveals  that 
the  ease  in  computation  offered  by  the  reduction  of  the  number  of  constraints  in  flatness-based 
methods  can  also  be  achieved  by  exploiting  spare  linear  algebra  in  solving  the  underlying 
optimization  problem  [39].  Further,  real-time  trajectory  generation  is  not  merely  driven  by  the 
number  of  constraints  or  optimization  variables  but  fundamentally  intertwined  with  optimization 
principles  such  as  convexity[40].  We  note  that  flatness-based  methods  belong  to  the  category  of 
direct  methods[l],  and  direct  comparison  between  flatness-based  methods  and  the  reduced-order 
formulation  is  not  feasible.  However,  with  the  given  scenario  of  terrain-masking,  moving  targets 
and  obstacles,  feasibility  study  on  the  structural  property  such  as  flatness  of  the  air  vehicle 
dynamics  will  further  shed  light  on  the  structure  of  the  air  vehicle  trajectory  optimization  and 
strengthen  the  understanding  on  the  viability  of  the  reduced-order  fonnulation. 
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Finally,  while  the  reduced-order  optimal  solution  in  a  certain  case  leads  to  a  trajectory 
that  can  not  be  tracked  with  the  full  aircraft  dynamics,  the  neural  network-based  autopilot  hedges 
the  unfollowable  portion  of  the  guidance  command  by  estimating  achievable  portion  of 
trajectories  by  aircraft  dynamics[41]  and  still  maintains  reasonable  flights.  For  example,  Figures 
5.25,  5.29,  and  5.32  shows  that  the  altitude  command  is  being  tracked  reasonably  well 
considering  significant  fluctuation  of  original  commands.  This  raises  a  question  if  guidance- 
command  hedging  mechanism  can  be  used  to  correct  the  feasibility  of  the  optimal  solutions  so 
that  a  resulting  flight  path  is  near  optimal  to  the  original  full-order  optimal  trajectory.  Some 
implications  provided  by  this  study,  such  as  simplification  of  re -planning  and  robustness  to  the 
modeling  error  that  will  destroy  the  optimality  anyway  if  the  aircraft  dynamics  are  uncertain, 
may  be  worth  further  study. 
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